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COMPUTAT  NAL  STUDIES  OF  LAMINAR  TO 
TI  LENCE  TRANSITION 


Introduction 

The  subject  of  laminar/turbulent  transition  is  of  fundamental  and  practical  importance  in  fluid 
mechanics.  An  indepth  knowledge  of  transition  mechanism  is  needed  not  only  for  boundary/shear 
layer  control  but  also  for  understanding  of  turbulence.  Depending  upon  the  state  of  the  boundary 
layer,  various  instability  mechanisms  such  as  Tollmien-Schlichting  (TS),  crossflow  and  Gdrtler  may 
be  operative.  The  present  day  transition  prediction  methodology  is  based  upon  linear  stability  theory 
in  the  form  of  method.  This  technique  works  well  in  a  specialized  parameter  space  involving 
extremely  “clean”  flows.  However,  depending  upon  the  flow  conditions,  various  instabilities  referred 
to  above  could  interact.  In  case  of  such  wave-interactions,  the  transition  prediction  methodology, 
based  solely  upon  linear  stability  theory,  would  fail.  Therefore,  it  is  of  interest  to  study  possible 
wave  interactions  and  devise  new  transition  prediction  criterion  under  such  circumstances. 

The  early  stages  of  breakdown  of  Tollmien-Schlichting  waves  to  turbulence  in  a  two-dimensional 
flat  plate  boundary  layer  are  now  relatively  well  understood  due  to  careful  laboratory  experiments, 
(KlebanofT  et  al.  1962,  Nishioka  et  al.  1980)  ingenious  numerical  simulations  (Zang  and  Hussaini 
1985,  Hussaini  1987,  Spalart  and  Yang  1987,  Fasel  et  al.  1987),  and  analytical  studies  (Herbert 
1988,  Kachanov  1987).  However,  our  knowledge  of  transition  mechanism  in  three-dimensional 
boundary  layers  and  flow  where  primary  instability  mechanism  is  other  than  TS  is  relatively  limited. 

A  swept  wing  boundary  layer  or  a  boundary  layer  on  a  body  at  an  angle  of  attack  is  s  -  ' 

both  crossflow  and  TS  type  instabilities.  While  the  former  results  due  to  inflectional  insta^''  ,  vf 
the  crossflow  velocity  profiles,  the  latter  is  a  viscous  instability  of  the  streamwise  profiles.  Another 
boundary-layer  where  two  types  of  instabilities  exist  is  that  formed  on  a  concavely  curved  plate. 
Here,  counter-rotating  steady  Gdrtler  vortices  form  due  to  centrifugal  instability.  Given  that  the 
Reynolds  number  is  high  enough,  TS  instability  may  also  be  present.  Nonlinear  development  of 
Gdrtler  vortices  and  Gdrtler/TS  interaction  is  again  a  problem  of  both  fundamental  and  potentially 
practical  importance. 

The  objective  of  the  present  study  is  to  understand  physical  mechanisms  leading  to  transition 
and  identify  any  new  instability  mechanisms  involved  in  these  flows.  A  numerical  study  is  carried 
out  to  investigate  how  a  steady  Gdrtler  vortex  breaks  down  to  turbulence.  The  study  reveals  various 
stages  involved  including  the  nonlinear  development  of  Gdrtler  vortex,  secondary  instabilities 
leading  to  the  development  of  unsteadiness  and  the  final  breakdown. 

In  order  to  investigate  crossflow  instability,  we  consider  a  model  swept-wing  boundary-layer  and 
study  the  development  of  stationary  as  well  as  traveling  disturbances  and  their  nonlinear  interac¬ 
tion.  It  is  also  found  that  prior  to  laminar-turbulent  transition,  the  three-dimensional  boundary 
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layer  is  subject  to  a  high-frequency  secondary  instability,  which  is  in  agreement  with  the  experi¬ 
ments. 

Two  approaches  have  been  used  in  the  above  studies.  For  nonlinear  evolution  and  wave-wave 
interaction,  we  utilize  parabolized  stability  equations  (PSE).  For  secondary  instability,  we  employ 
two-dimensional  (2D)  eigenvalue  approach. 

The  final  report  is  divided  into  three  parts; 

A.  On  the  Breakdown  of  Gortler  Vortices:  Nonlinear  Development  and  Secondary  Instabilities 
This  describes  nonlinear  evolution  of  Gortler  vortices  and  their  breakdown  via  various  sec¬ 
ondary  instability  mechanisms. 

B.  Crossflow  Disturbances  in  Three-Dimensional  Boundary  Lavers:  Nonlinear  Development.  Wave 
Interaction  and  Secondary  Instability 

This  describes  nonlinear  evolution  of  crossflow  disturbances,  wave-wave  interaction  and  a  high- 
frequency  secondary  instability  prior  to  laminar-turbulent  breakdown. 

C.  On  the  Nature  of  PSE  Approximation 

Since  the  PSE  approach  is  used  in  the  above  physical  problems,  it  is  important  to  investigate 
the  mathematical  nature  of  the  approximation.  In  this  part,  we  explore  the  parametric  boundaries 
within  which  PSE  approximation  is  valid  and  results  in  a  stable  numerical  solution.  This  knowledge 
is  used  when  we  employ  this  approach  in  the  above  problems. 
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PART  A. 

On  the  Breakdown  of  Gortler  Vortices:  Nonlinear  Development 

and  Secondary  Instabilities 
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Abstract 

Nonlinear  development  of  stationary  Gdrtler  vortices  in  an  incompressible  boundary  layer  is 
studied  by  solving  the  parabolic  partial  differential  equations.  It  is  found  that  due  to  the  pumping 
action  of  the  steady  counter-rotating  vortices,  wall  shear  stress  decreases  rapidly  in  the  peak  plane 
but,  unlike  the  result  of  Hall  (1990),  it  does  not  go  to  zero  and  there  is  no  flow  reversal.  Instead,  as 
the  (jdrtler  vortex  saturates  the  mean  flow  correction  mode  dominates  and  causes  the  wall  shear  in 
the  peak  plane  to  begin  to  rise  again.  A  highly  distorted  mean  flow  field  is  set  up  due  to  the  vortex 
action  where  the  streamwise  velocity  u  depends  strongly  not  only  on  y  (wall  normal)  but  also  on  the 
z  (spanwise)  coordinate.  The  inviscid  instability  of  this  flow  field  is  analyzed  by  solving  the  two- 
dimensional  eigenvalue  problem  associated  with  the  governing  partial  differential  equation.  It  is 
found  that  the  flow  field  is  subject  to  odd  and  even  (with  respect  to  the  (Jdrtler  vortex)  unstable 
modes.  The  odd  mode  which  was  also  found  by  Hall  &  Horseman  (1991)  is  initially  more  unstable. 
However,  there  exists  an  even  mode  which  has  higher  growth  rate  further  downstream.  The 
nonlinear  development  of  these  secondary  instability  modes  is  studied  by  solving  the  (viscous)  partial 
differential  equations  under  a  parabolizing  approximation.  The  odd  mode  leads  to  the  well-known 
sinuous  mode  of  breakdown  while  the  even  mode  leads  to  the  horse-shoe  type  vortex  structure.  This 
helps  explain  experimental  observations  that  Gdrtler  vortices  breakdown  sometimes  by  sinuous 
motion  and  sometime  by  developing  a  horse-shoe  vortex  structure.  The  details  of  these  breakdown 
mechanisms  are  presented. 

1.  Introduction 

Two-dimensional  boundary-layer  flow  over  a  concavely  curved  wall  is  subject  to  Gdrtler 
instability  due  to  the  action  of  centrifugal  force  and  results  in  the  formation  of  counter-rotating 
streamwise  vortices.  (Jdrtler  vortices  play  a  dominant  role  in  boundary-layer  transition  in  many 
aerodynamic  flows  such  as  on  turbine  blades  and  supersonic  nozzle  walls  (e.g.,  Beckwith  et  al.  1984). 
Due  to  their  technological  importance,  Gdrtler  vortices  have  been  the  subject  of  a  number  of 
investigations  (for  a  recent  review,  see  Floryan,  1991).  (Jdrtler  vortices  are  steady  and  the  question 
how  they  might  breakdown  to  turbulent  motion  is  a  problem  of  fundamental  interest  in  fluid 
mechanics.  This  problem  may  also  serve  as  a  model  for  the  longitudinal  vortices  in  turbulent 
boundary  layers.  In  this  work,  we  will  study  the  nonlinear  development  of  (Jdrtler  vortices,  their 
linear  secondary  instability  characteristics  and  the  nonlinear  growth  of  two  important  modes  of 
secondary  instability  up  to  the  breakdown  stage. 

Experimental  investigations  have  revealed  two  distinct  types  of  secondary  instabilities  when  the 
primary  instability  (the  (Jdrtler  vortex  )  is  sufficiently  developed.  Bippes  (1978)  made  detailed 
observations  of  the  Gdrtler  vortex  breakdown  using  the  hydrogen-bubble  visualization  technique  in 
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the  Gortler  number  (Gg)  range  of  approximately  3  to  9  (based  upon  momentum  thickness  d)  where 
Gg  is  defined  as 

(1.1) 

being  the  surface  curvature  and  Rg  the  Reynolds  number  based  upon  momentum  thickness. 
Bippes  found  that  the  initial  amplification  of  the  (jdrtler  vortices  agreed  with  linear  theory  and, 
later,  sinuous  oscillations  developed,  which  ultimately  led  to  turbulence.  Aihara  &  Koyama  (1981) 
conducted  flow  visualization  studies  as  well  as  hot-wire  measurements  of  Gortler  vortices  in  Gortler 
number  (Gg)  range  between  7.7  and  15.  They  found  that  a  different  type  of  secondary  instability, 
i.e.,  the  horse-shoe  vortex  t}rpe  (also  called  the  varicose  instability),  was  responsible  for  transition. 
Ito  (1985)  also  found  this  ssnnmetric  mode  of  breakdown  in  his  experiment  conducted  in  the  Gortler 
number  range  between  5.5  and  12.4.  Swearingen  &  Blackwelder  (1987)  (referred  to  as  SB  hereafter) 
studied  the  Gdrtler  vortices  using  smoke-wire  and  hot-wire  techniques  in  the  (Jdrtler  number  range 
between  0.5  and  10.  They  observed  that  both  the  sinuous  and  the  horse-shoe  types  of  secondary 
instabilities  were  present  in  the  transition  process  and  found  that  the  sinuous  mode  was  the  stronger 
of  the  two.  In  their  experiment,  the  unsteady  secondary  instability  fluctuations  correlated  better 
with  the  spanwise  velocity  gradients  than  with  the  normal  velocity  gradient.  Unsteady  motion  in 
Gortler  vortices  was  also  observed  by  Peerhossaini  &  Wesfried  (1988). 

Numerical  simulations  were  carried  out  by  a  number  of  researchers  for  Gortler  vortices  under 
the  same  conditions  as  the  experiment  of  SB.  Sabry  &  Liu  (1991)  studied  the  key  features  of  steady 
Gortler  vortex  development  by  using  a  temporal  analogy  in  which  the  growth  of  time-dependent 
streamwise  vortices  are  related  to  the  spatial  case  through  a  chosen  advection  velocity.  Good 
agreements  with  SB  were  found  before  the  unsteady  oscillations  became  important.  Liu  & 
Domaradzki  (1993)  used  a  similar  approach  and  analyzed  the  steady  and  unsteady  nonlinear 
evolution  of  Gortler  vortices.  Temporal  Navier-Stokes  approach  was  earlier  used  by  Malik  (1986) 
and  Malik  &  Hussaini  (1990)  to  study  (Jdrtler/ToUmien-Schlichting  wave  interaction.  Since  the 
physical  problem  is  spatial,  temporal  approach  can,  at  best,  provide  qualitative  results  and  there  is 
no  a  priori  justification  for  this  approximation.  Spatial  simulations  of  the  (Jortler  vortex  were 
performed  by  Hall  (1988,  1990)  and  by  Hall  &  Horseman  (1991).  Agreement  with  the  experiment  of 
SB  was  obtained  for  the  early  stages  of  (Jdrtler  vortex  development.  However,  Hall  (1990)  and  Hall 
&  Horseman  (1991)  found  that  the  wall-shear  in  the  "peak-plane"  crossed  zero  to  become  negative  at 
some  downstream  location,  indicating  the  existence  of  a  region  of  reverse  flow.  They  attributed  the 
rise  in  wall-shear  in  the  experiment  of  SB  to  the  nonlinear  interaction  of  unsteady  oscillations  and 
the  steady  Gortler  vortex.  However,  spatial  calculations  carried  out  by  Lee  &  Liu  ( 1992)  did  not  find 
any  flow  reversal. 


5 


Various  theoretical  attempts  at  the  secondary  instability  of  the  Gdrtler  vortex  illuminated  the 
important  role  the  sinuous  and  the  varicose  modes  play  in  transition  to  turbulence.  Sabry,  Yu  &  Liu 
(1990)  used  1-dimensional  local  inflectional  velocity  profiles  to  analyze  the  secondary  instability 
characteristics  of  Gdrtler  vortices  and  found  that  the  sinuous  type  disturbance  would  prevail  over  the 
varicose  type.  Hall  &  Horseman  (1991)  derived  a  partial  differential  equation  governing  the  inviscid 
secondary  instability  for  a  mean  flow  which  varied  strongly  in  two  directions.  They  identified  the  2- 
dimensional  odd  and  even  eigenfunctions  of  this  equation  as  representing  the  sinuous  and  varicose 
types  of  secondary  instability  of  Grdrtler  vortex,  respectively.  They  also  found  that  the  odd  mode 
grows  faster  than  the  even  mode.  In  the  numerical  simulation  of  Liu  &  Domaradzki  (1993),  the 
sinuous  mode  was  also  found  to  be  stronger  than  the  varicose  mode.  Secondary  instability  analysis 
of  Gdrtler  and  crossflow  disturbances  was  made  by  Malik  and  Li  (1993a). 

All  these  numerical  investigations  of  secondary  instability  seem  to  point  to  one  fact,  i.e.,  the 
sinuous  mode  is  the  dominant  mode  and  is  chiefly  responsible  for  the  transition  to  turbulence.  Why 
do  some  experiments  (  e.g.,  Aihara  &  Koyama,  1981)  show  the  presence  of  only  the  varicose  mode? 
Examining  the  above  cited  numerical  results  closely,  we  find  that  most  of  these  calculations  were 
carried  out  for  either  a  limited  range  of  wavelengths  or  a  limited  number  of  streamwise  locations.  In 
Hall  &  Horseman  (1991),  for  example,  the  secondary  instability  calculations  were  computed  almost 
exclusively  at  x  =  100  cm.  In  Liu  &  Domaradzki  (1993),  the  computational  box  had  streamwise 
dimension  of  either  2  cm  or  2.2  cm,  which  essentially  fixed  the  streamwise  wavelength. 
Furthermore,  a  comparison  of  the  "mushroom"  structure  obtained  by  Hall  &  Horseman  (1991)  on  the 
one  hand  and  Lee  &  Liu  (1992)  and  Liu  &  Domaradzki  (1993)  on  the  other  shows  that  the  basic  flow 
state  used  for  secondary  instability  computations  are  different.  In  the  results  of  Hall  &  Horseman 
(1991),  the  "mushroom"  lacks  the  very  thin  "stem”  shown  in  the  results  of  other  researchers. 

In  this  paper,  we  will  first  study  the  nonlinear  development  of  the  (Jdrtler  vortex  and  analyze 
various  modes  of  secondary  instability.  Rayleigh's  criterion  for  the  inviscid  instability  of  one¬ 
dimensional  velocity  profile  is  well-known.  However,  its  extension  to  two-dimensional  flow  has  not 
been  considered.  We  will,  therefore,  derive  a  necessary  condition  for  inviscid  instability  of  mean 
flows  which  strongly  depend  upon  two  space  variables.  We  also  study  nonlinear  evolution  of  the 
secondary  instability  modes  using  parabolized  stability  equations  (PSE).  Section  2  describes  the 
parabolized  stability  equation  formulation  which  describes  all  stages  of  the  vortex  development. 
Section  3  deals  with  the  nonlinear  development  of  steady  Gdrtler  vortices  while  sections  4  &  5 
anal3rzes  the  linear  secondary  instability  and  its  nonlinear  evolution,  respectively. 

2.  Problem  Formulation  for  Nonlinear  Steady  and  Unsteady  Disturbances 

We  consider  a  two-dimensional  zero  pressure  gradient  boundary-layer  flow  over  a  concave 
surface  whose  constant  radius  of  curvature  =  The  streamwise,  wall-normal  and  the 
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gpanwise  coordinates  are  denoted  as  x  =  X/^o,  y-Y/f^  and  z-ZI(q,  respectively  (y  =  0  denotes 
the  wall),  where  the  length  scale  will  be  prescribed  later. 

Let  the  x,  y  and  z  components  of  the  velocity  and  pressure  be  given  by 


=  f/,  y)  +  u(x, y, /),  V(jc, y)  +  v{x, y, z, t),  0  +  w{x,  y,z  t)] 

=  pU'^{P  +  p{x,y,z,t)) 

where  the  superscript  t  represents  a  dimensional  quantity  and  is  the  velocity  scale.  Here  U  and 
V  are  mean  flow  velocity  components  obtained  by  solving  the  Blasius  equation  whereas  u,  v,  w 
represent  the  perturbation  velocity  components  in  x,  y,  z  directions,  respectively.  Similarly,  P  and  p 
represent  the  mean  and  perturbation  pressures.  We  assume  that  Reynolds  number,  R,  is  large  and 
that  the  radius  of  curvature  is  much  larger  than  the  boundaiy-layer  thickness,  S  (i.e.,  k^S  «  1).  In 
this  case,  if  v  =  K^io,  the  equations  governing  the  perturbation  quantities  are 


du  dU  .,du  dU  ..  K  dp  1  „2  »r 

—  +  U  —  +  u—-  +  V  —  +  v—-  +  K[Uv  +  Vu)  +  -f--  —  V‘'u  =  N, 
dt  dx  dx  dy  dy  ^  ’  dx  R  ^ 

du  ..du  dV  ^^du  dV  dp  1  „2 

dt  dx  dx  dy  dy  dy  R 


du) 

dt 


+  f/^  +  V^  +  ^-lv2u,  =  iV. 


dx 


dy  dz  R 


du  d)  dw  „ 

where  and  represent  the  nonlinear  terms 

du  du  du 

Ni=-u—-v—-w~-kuu 
dx  dy  dz 


du  du 


d) 


•  T  W  VW  WV’ 

No  =-U—-U—-W—  +  KU 


dx  dy 


dw  du 


dz 


dw 

''“'-“a-"*-"’* 


and 


d^  d^  d^  d 

d)r  dy^  dz^  dy 


The  boundary  conditions  are 

u  =  v  =  w  =  0  at  y  =  0  and  u-»0,  i;->V„(x),  ie-)0,  as  y~^°° 


(2.1) 

(2.2) 

(2.3) 

(2.4) 
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where  V„  signifies  a  nonzero  value.  Periodic  boundary  conditions  are  imposed  in  the  spanwise 
direction.  Here,  (q  and  Lf,  are  constants  so  as  to  define  Reynolds  number  R  as 


V 


vXa 


where  the  length  scale  (q  =  being  the  location  (dimensional)  of  a  reference  streamwise 

station,  and  v  the  kinematic  viscosity.  Another  important  parameter  which  is  a  measure  of  the  wall 
curvature  is  the  Gdrtler  number  G  defined  as  G  =  .  The  reason  why  v  does  not  go  to  zero 

outside  the  boundary  layer  is  that  the  vertical  velocity  vanishes  for  all  Fourier  modes  (see  (2.5) 
below)  except  the  mean-flow  correction  mode. 

We  use  the  method  of  parabolized  stability  equations  (PSE)  for  our  computations.  Here,  we 
follow  the  formulation  of  Malik  &  Li  (1993b)  (see  also  Malik  et  al.  (1994))  and  let  0  =  {u,v,w,p)  be  the 
disturbance  vector  and  assume  that  the  disturbance  takes  the  form 


mn  (x,  y)  exp]  1 1"  a„„  (|)cf  =  +  inpz  -  im  ait 

m  n  I  •'*0  J 


(2.5) 


where  a^n  and  p  are  the  x  and  z  wave  numbers,  o)  is  the  perturbation  frequency  and  is  the 

amplitude  function  for  the  mode  (ma),n/J).  Substituting  Eq  (2.5)  into  Eqs  (2.1)  to  (2.4),  we  obtain  a 

set  of  equations  with  0„,„  and  a„„  as  unknowns.  Since,  there  are  now  more  unknowns  (namely, 

a„„)  than  the  equations,  another  condition  is  needed  for  the  closure  of  the  system.  Since  the  basic 

flow  is  slow-varying  in  the  streamwise  direction,  a  condition  on  is  imposed  such  that  most  of  the 

waviness  and  growth  of  the  perturbation  is  absorbed  into  the  exponential  function  in  Eq  (2.5), 

making  the  amplitude  function  0„,„  slowly-varying  with  respect  to  x  .  The  terms  containing  — 

dx 

can  thus  be  dropped  and  the  only  second  derivatives  left  in  the  governing  equations  are  those  with 
respect  to  y.  These  new  stability  equations  are  parabolized  in  the  sense  of  parabolized  Navier-Stokes 
(PNS)  equation  for  mean  flow  computations.  The  condition  for  choosing  a„„  and  minimizing  the 
streamwise  variation  of  the  amplitude  function  can  take  several  forms.  In  the  present  work,  we 
choose  a„„  to  be  such  that  the  following  integral  vanishes. 


/a»  «•  d 

\u  ,v  ,w  )-^ 

V 

^  'dx 

dy  =  0 


(2.6) 


where  *  denotes  complex  conjugate.  The  PSE  can  be  written  in  matrix  form  as 


(2.7) 
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(2.8) 


where  the  coefficient  matrices  contain  the  Blasius  flow  quantities  as  well  as  a^„,  P  and  ax  Eq  (2.8)  is 
a  general  form  of  Eq  (2.6).  The  matrix  operators  L^-L^  and  are  given  in  Appendix  I.  The 
boundary  conditions  are 

“mn  =  =  0  at  y  =  0  (2.9a) 

“mn. (except  for  m  =  n  =  0), UI„„ -^0  as  y->°°  (2.9b) 

No  boundary  condition  is  requried  for  v{°°)  when  m  =  n  =  0. 

We  discretize  the  PSE  using  discrete  Fourier  transforms  in  the  spanwise  direction  and  in  time. 
In  the  direction  normal  to  the  wall,  we  use  the  fourth-order  compact  difference  scheme  (Malik, 
Chuang  &  Hussaini,  1982)  which  requires  that  (2.7)  be  written  as  a  system  of  first-order  equations. 
Numerical  computation  starts  at  some  streamwise  location  Xq  where  velocity  components  are 
prescribed  for  a  given  wavenumber  a„„  ;  the  velocities  and  pressure  at  Xq  calculated  using 

backward  Euler  discretization.  If  Eq  (2.8)  were  not  satisfied,  a  new  would  be  chosen  and  the 
equations  solved  again.  This  iterative  process  continues  until  Eq  (2.8)  is  satisfied,  and  the 
computation  proceeds  to  the  next  streamwise  location.  During  this  iterative  process,  nonlinear  terms 
are  also  updated  and  one  makes  sure  that  they  are  converged  before  the  solution  proceeds 
downstream. 

3.  Nonlinear  Development  of  Steady  Gortler  Vortices 

In  this  section,  we  analyze  the  nonlinear  development  of  streamwise  stationary  Gdrtler  vortices 
(a)  =  a  =  0).  In  the  limit  /2  ->  <» ,  and  x  0  with  G  held  fixed,  and  by  rescaling  the  dependent  and 
independent  variables  (V  =  0(1/ It)U,(v,u;)  =  0(1/  R)u,y  =  0(1/ Ii)x)  the  parabolic  equations  derived 
by  Hall  (1983,  1988)  can  be  recovered  from  (2. 1-2.4).  In  Hall  (1988),  a  further  step  is  taken  to 
eliminate  the  spanwise  velocity  and  the  pressure  from  the  linear  terms  in  Eq  (2.7),  resulting  in  a 
coupled  system  of  fourth-order  and  second-order  equations.  However,  we  will  solve  Eq  (2.7)  directly 
in  the  primitive  (u,v,w,p)  formulation,  except  that  the  condition  imposed  on  a  (2.8)  is  not  applied 
since  a  is  identically  0. 

The  flow  parameters  used  in  the  present  analysis  are  taken  from  the  experiment  of  SB.  The 
radius  of  curvature  of  the  concave  surface  is  320  cm  and  free  stream  velocity  is  500  cm/s.  The 
streamwise  range  of  interest  in  our  analysis  lies  approximately  in  10  <  X  (cm)  <  120,  in  which  the 
Reynolds  number  based  on  the  distance  from  the  leading  edge  (Re  =  i?^)  ranges  from  3.3x10^  to 
4x10®  and  the  Gortler  number  Gg  ranges  from  1.3  to  8.3.  The  Blasius  equations  are  solved  to 
obtain  the  basic  flow  {U,V)  for  the  vortex  stability  analysis.  The  calculation  is  started  at  Xq  =10  cm 
for  the  disturbance  wavelength  =  1.8  cm  and  initial  amplitude  of  i2  =  0.0187  C/^  estimated  from 
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the  experimental  data  of  SB.  At  this  stage,  our  aim  is  to  analyze  the  steady  nonlinear  Gdrtler 
vortices;  therefore,  we  only  allow  zero-frequency  modes  to  be  present.  Unsteady  modes  will  be  taken 
into  account  later  when  we  compute  the  development  of  secondary  instabilities.  The  number  of 
spanwise  Fourier  modes  used  in  the  z-direction  is  11  (i.e.,  n  goes  from  -10  to  10  in  (2.5)),  and  the 
streamwise  marching  step-size  is  0.82158  cm.  Therefore,  approximately  130  marching  steps  are 
taken  in  the  streamwise  direction.  The  number  of  wall-normal  steps  is  121.  The  solution  was  tested 
by  changing  the  number  of  grid  points  and  it  was  found  to  be  grid-independent. 

Attempt  is  made  to  compare  with  the  experiment  of  SB  in  this  paper.  A  few  words  of 
explanation  are  needed  in  order  to  clarify  the  manner  in  which  comparisons  are  made.  There  is 
uncertainty  in  the  experiment  of  SB  (and  in  any  experiment,  in  general)  with  regards  to  the 
spanwise  wavelength  of  the  Gdrtler  vortices.  As  noted  by  SB,  this  wavelength  varies  across  the  span 
and  its  statistical  average  is  2.3  cm.  Th  detailed  measurements,  however,  are  made  for  a  pair  of 
vortices  with  a  spanwise  wavelength  of  about  1.8  cm  which,  in  fact,  is  close  to  the  most  amplified 
Gdrtler  vortex  according  to  linear  theory.  There  is  also  a  fa  r  amount  of  scatter  in  the  data  in  the 
"peak"  region  of  the  vortices  for  localized  quantities  such  as  the  wall-shear.  Our  computation  is 
performed  with  a  constant  spanwise  wavelength,  which  is  fixed  at  1.8  cm.  We  could  choose  different 
wavelengths  and  initial  amplitudes  to  obtain  a  whole  range  of  results  from  which  the  best 
comparison  with  experiment  can  be  found,  but  we  would  have  little  to  gain  from  this.  As  Hall  & 
Horseman  (1991)  pointed  out  that  the  inherent  non-uniqueness  of  the  Gdrtler  problem  might  be 
present  in  the  experiment  as  well  and  the  exact  features  of  the  experiment  itself  might  not  be 
precisely  reproducible.  Consequently,  our  comparisons  with  the  experiment  are  confined  to  the 
qualitative  features  and  trends  of  the  Gdrtler  vortex  development  which  are  relatively  insensitive  to 
moderate  variations  in  parameters  such  as  wavelengths  and  initial  amplitudes. 

The  experiment  of  SB  (as  well  as  others)  produced  "mushroom-like"  structures  for  the 
streamwise  velocity  due  to  the  pumping  action  of  the  counter-rotating  vortices.  The  contours  of  the 
streamwise  velocity  at  various  downstream  locations  computed  in  the  present  analysis  are  shown  in 
Figure  1.  In  the  early  stages  of  the  development,  the  amplitude  of  u  perturbation  is  small  and  the 
velocity  contours  show  a  wavy  spanwise  structure.  As  the  Gdrtler  vortices  gather  strength  at 
relatively  large  distances  downstream  the  same  "mushroom"  structures  observed  in  the  experiment 
of  SB  are  clearly  seen.  The  regions  in  the  neighborhood  of  the  centerlines  of  the  "mushrooms"  are 
referred  to  as  "peak"  regions  where  the  streamwise  velocity  is  relatively  low;  and  the  regions  between 
the  "mushrooms"  are  referred  to  as  "valley"  regions  where  the  streamwise  velocity  is  relatively  high. 
The  streamwise  velocity  profiles  at  the  peak  and  the  valley  are  shown  in  Figure  2.  It  will  be  shown 
later  that  the  high  shear  layer  region  up  in  the  peak  plane  will  become  subject  to  a  particular  mode 
of  secondary  instability.  The  spanwise  variation  of  the  streamwise  velocity  component  at  fixed  y  is 
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given  in  Figure  3.  Again,  the  inflected  profiles  will  be  subject  to  another  mode  of  secondary 
instability. 

The  energy  in  each  Fourier  mode  is  shown  in  Figure  4.  Here,  energy  is  defined  as 


E„  = 


(l“nr  ^\l^nf)dy 


where  n>0,  and 


(3.1a) 


(3.1b) 


Initially,  the  (jdrtler  vortex  (mode  1)  develops  approximately  linearly,  but  later  on  it  begins  to 
saturate.  Due  to  nonlinearity,  higher  harmonics  (mode  2,3...)  and  mean  flow  distortion  mode  (mode 
0)  are  generated.  Among  these,  first  the  mode  0  and  2  have  about  the  same  energy,  but  at 
approximately  X  =  75  cm,  the  energy  in  the  mean  flow  correction  mode  takes  over  that  of  the 
fundamental,  and  hence,  the  mean  flow  correction  mode  becomes  the  dominant  mode,  as  also  found 
by  Hall  and  Lakin  (1988)  for  the  small  wavelength  Gdrlter  vortices.  This  has  important  physical 
consequences,  which  will  be  explained  later.  Farther  downstream,  the  energy  in  each  mode  begins  to 
level  off,  approaching  saturation.  Although  the  energy  in  modes  2  and  higher  remains  much  lower 
than  that  of  the  fundamental,  they  all  are  important  in  the  development  of  the  narrow  neck 
structure  of  the  mushroom  seen  in  Figure  1(d). 

We  now  consider  the  wall  shear  in  the  peak  and  the  valley  plane.  FigTure  5  shows  the  computed 
normal  wall  shear  in  these  regions,  together  with  the  experimental  data  of  SB.  Initially,  the  shear  at 
the  peak  and  at  the  valley  decreases  and  increases,  respectively,  due  to  the  slow-down  and  speed-up 
of  streamwise  velocity  at  the  respective  locations.  Farther  downstream,  the  shear  at  the  peak  turns 
around  and  starts  to  increase  at  approximately  X  =  70.  It  has  been  argued  by  Hall  (1990)  that  the 
rise  in  shear  is  solely  due  to  nonlinear  interaction  of  steady  and  unsteady  disturbances,  without 
which  the  shear  would  continue  to  decrease  and  eventually  leads  to  flow  reversal  near  the  wall. 
However,  in  SB,  the  amplitude  of  the  unsteady  fluctuations  at  this  stage  is,  at  least,  one  order  of 
magnitude  smaller  than  that  of  the  steady  (jdrtler  vortex  and  it  is  unlikely  that  such  drastic  change 
in  flow  character  is  caused  by  these  unsteady  fluctuations.  The  normal  wall  shear  at  the  peak  is 
larger  compared  with  the  experiment  of  SB.  Using  a  spanwise  wavelength  of  2.3  cm  instead  of  1.8 
cm  gives  a  much  better  comparison.  We  point  out  that  Lee  &  Liu  (1992)  used  a  lower  initial 
amplitude  for  the  wavelength  for  =  1.8  cm  and  showed  better  agreement  with  experiment  for  wall 
shear.  We  point  out  that  in  our  calculations,  the  increase  in  the  shear  stress  is  independent  of  the 
initial  conditions  and  the  location  where  they  are  imposed. 
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Now  we  explore  the  reason  why  the  wall  shear  at  the  peak  region  turns  around  and  starts  to 
increase  at  approximately  X  =  70  cm  for  A,  =  1.8  cm.  This  location  is  very  close  to  the  point  where 
the  energy  in  the  mean  flow  correction  mode  takes  over  that  of  the  fundamental  (X  =  72  cm).  Also 
plotted  in  Figure  5  is  the  streamwise  mean  wall-shear,  which  is  entirely  due  to  the  mean  flow 
correction  to  which  all  the  harmonics  contribute.  The  mean  wall  shear  at  large  distances 
downstream  has  rather  large  positive  amplitude.  Therefore,  the  explanation  for  the  wall  shear  turn¬ 
around  lies  in  the  fact  that,  while  the  Gdrtler  vortices  try  to  slow  down  the  fluid  in  the  peak  region, 
the  mean  flow  correction  tries  to  accelerate  it.  As  the  mean  flow  correction  becomes  more  and  more 
significant  and  eventually  dominant,  the  wall  shear  has  a  large  contribution  from  the  mean  flow 
correction  mode,  which  overpowers  the  other  modes  to  produce  an  increase  in  the  streamwise 
velocity  in  the  peak  plane,  resulting  in  the  increase  in  wall  shear.  Figure  6  shows  the  computed 
streamwise  velocity  profiles  along  the  peak  region  at  various  streamwise  locations.  It  clearly  shows 
(see  also  Figure  2)  that  even  though  the  profiles  are  highly  inflectional,  the  near  wall  behavior  of  the 
velocity  would  not  give  rise  to  negative  shear.  The  domination  of  the  mean  flow  correction  mode  also 
explains  why  the  shear  stress  cannot  go  to  zero  and  cause  flow  reversal. 

The  evolution  of  displacement  thickness  is  shown  in  Figure  7.  The  computed  results  agree  well 
with  the  experiment  of  SB  up  to  about  X  =  100  cm  where  unsteady  oscillations  in  the  experiment 
become  important.  The  displacement  thickness  decreases  at  the  valley  because  of  the  high  speed 
fluid  brought  down  to  the  region  near  the  wall  by  the  Gdrtler  vortices.  In  the  peak  region,  the 
displacement  thickness  increases  much  faster  than  the  corresponding  displacement  thickness  of  the 
Blasius  boundary  layer  because  of  the  low  speed  fluid  brought  into  the  upper  regions  of  the  boundary 
layer. 

The  present  computations  showed  many  of  the  qualitative  features  of  the  experimentally 
observed  Gdrtler  vortices.  However,  for  an  understanding  of  the  breakdown  of  these  vortices  due  to 
secondary  instabilities,  we  have  to  consider  unsteady  modes.  We  will  use  the  steady  Gdrtler  vortex 
flow  field  generated  above  as  the  basic  flow  state  for  our  secondary  instability  ansdysis  given  in  the 
next  section.  Thus,  if  u^,  Vq,  Wq  are  the  perturbation  velocities  due  to  the  Gdrtler  vortex,  the  new 
mean  flow  whose  stability  will  be  analyzed  in  the  next  section  is 


U=U  +Uq 

(3.2a) 

U  =  V  +  Vq 

(3.2b) 

w  =  W  +  wo 

(3.2c) 
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4.1 


We  assume  that  the  stream  wise  variation  of  the  new  mean  flow  (u,U,w)  is  small  compared  with 
the  wavelength  of  the  secondary  disturbance.  This  assumption  can  be  justified  a  posteriori  from  the 
results.  Therefore,  the  linear  secondary  oscillation  can  be  written  as 

=  (4.1) 

where  ^  =  {u,,v„w„Pg)  and  \ff  =  {u„v,,Wg,Pgy,  a  is  the  wavenumber,  o)  is  the  oscillation  frequency, 
and  the  subscript  s  indicates  secondary  instability.  Here,  we  consider  only  the  temporal  instability; 
therefore  a  is  real  and  (o  is  complex.  The  flow  field  is  unstable  to  disturbances  if  tUj  >  0. 

We  note  that  v(y,z)  and  w{y;i)  are  much  smaller  than  u(y,z).  Then,  following  Hall  & 
Horseman  (1991),  the  equations  governing  the  linear  secondary  instability  are 

.  .  30,  dw,  _  , . 

tow, +— ^  +  — ^  =  0  (4.2) 

3y  dz 

ia{u  -  c)t2,  =  -top.  (4.3) 

ia{u-c)v,~-^  (4.4) 

cfy 

ia{u~c)w,=-^^  (4.5) 

az 

where  c  =  (o/a  is  the  (complex)  disturbance  phase  velocity.  Eliminating  a,,  o,,  ta,,  we  obtain  the 
equation  governing  the  secondary  pressure  oscillation  (after  dropping  subscript  s): 

(4.6, 

ok  J  u-c  u-c 


The  boundary  conditions  are 
y  =  0,  p^(y,z)  =  0 

p(y,2)-»0 


(4.7a) 


(4.7b) 


p{y,2)  =  p{y,z  +  x^)  (4.8) 

where  A,  is  the  (jdrtler  vortex  wavelength. 

Equations  (4.6-4.8)  constitute  an  eigenvalue  problem  which  is  solved  by  using  a  Chebyshev 
collocation  method  in  the  y  direction  and  a  Fourier  collocation  method  in  the  z  direction  with 


appropriate  grid  stretchings  in  both  directions  to  concentrate  more  collocation  points  in  regions  of 
high  gradients.  Furthermore,  since  the  basic  flow  state  is  symmetric,  the  eigenfunctions  can  be  split 
into  families  of  even  and  odd  modes.  For  the  even  mode,  i^y,z)  =  p{y,-z)  and  for  the  odd  mode, 
p{y,z)  =  -p{y,-z).  Taking  advantage  of  the  sjnmmetry  conditions,  we  reduce  the  size  of  the  resulting 
discretized  system  by  approximately  half.  The  discretized  system  can  be  represented  in  the  form 

A<p  =  aB(p  (4.9) 

where  is  a  diagonal  matrix  and  A  is  a  square  matrix  of  size  Ny{N^I2  + 1),  where  Ny  and  are 
the  number  of  collocation  points  in  y  and  z  directions,  respectively.  This  eigenvalue  problem  is 
solved  by  the  QR  method  which  3rields  all  the  eigenvalues  of  the  discretized  system  (4.9). 
Throughout  the  computations,  we  use  Ny  =  85  and  N^  =  32. 

4.2  Results  of  Secondary  Instability  Calculations 

During  the  early  stage  of  the  Gfirtler  vortex  development,  the  quantity  is  negative 
everywhere  in  the  (y,z )-plane.  At  approximately  X  =  40  cm,  a  region  of  positive  begins  to  appear 
near  the  wall.  Therefore,  according  to  the  necessary  condition  for  instability  (  i.e.  V^u  =  0 
somewhere  in  the  flow  field),  the  (Sdrtler  vortex  for  the  initial  conditions  prescribed  in  §3  is  stable  to 
secondary  disturbances  for  X  <  40  cm.  We  note  that  we  have  not  shown  whether  the  above  condition 
is  also  a  sufficient  condition  for  instability. 

We  start  our  2-D  eigenvalue  computations  at  X  =  65  cm,  where  the  high  frequency  oscillations 
are  moderately  unstable.  Hence,  we  avoid  the  difficulty  associated  with  the  singularity  due  to 
neutral  disturbances.  The  growth  rate  variation  with  the  streamwise  wavenumber  at  various 
streamwise  locations  between  X  =  65  cm  and  X  =100  cm,  normalized  with  scales  at  X  =  10  cm  ,  is 
shown  in  Figure  8.  The  general  trend  is  that  the  secondary  oscillations  become  more  unstable  as  the 
Gldrtler  vortices  become  stronger  downstream.  The  maximum  growth  rate  at  each  streamwise 
location  occurs  at  streamwise  wavenumbers  approximately  between  0.2  and  0.3,  corresponding  to 
wavelengths  between  1.2  and  1.7  cm.  The  Blasius  boundary  layer  thickness  in  the  absence  of  the 
Gdrtler  vortices  in  the  range  between  X  =  65  cm  and  X  =  100  cm  is  approximately  between  0.7  and 
0.9  cm.  This  shows  that  the  wavelength  of  the  secondary  instabilities  is  of  the  order  of  boundary- 
layer  thickness.  Therefore,  our  assumption  that  the  basic  flow  state  variation  is  negligible  over  the 
distance  of  one  wavelength  is,  indeed,  justified.  We  can  visually  extrapolate  the  growth  rate  curves 
and  see  that,  in  this  streamwise  range,  the  highest  wave  number  where  secondary  instability  occurs 
is  approximately  0.5,  corresponding  to  a  wavelength  of  0.69  cm. 

We  now  consider  the  variation  of  the  maximum  growth  rate  of  the  secondary  instabilities  with 
streamwise  distance.  In  order  to  show  that  the  secondary  instability  grows  much  faster  than  the 
(jTortler  vortex,  we  convert  the  temporal  growth  rate  to  the  spatial  one  by  using  group  velocity 
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transformation.  The  use  of  group  velocity  to  transform  temporal  growth  to  spatial  growth  is  well 
known  for  boundary-layer  instabilities  (Gaster  1962,  Nayfeh  &  Padhye  1979).  The  value  of  group 
velocity  dm^  / da  varies  with  x  and  lies  in  the  range  of  0.6  to  0.72.  Figure  9  shows  the  spatial 
maximum  growth  rate  of  the  most  unstable  even  and  odd  modes,  as  well  as  the  spatial  growth  rate  of 
the  Gdrtler  vortex.  The  odd  mode  begins  to  become  significant  from  approximately  X  =  65  cm,  and 
the  even  mode  roughly  from  X  =  75  cm.  An  important  feature  we  discover  here  is  that,  although  the 
odd  mode  is  the  first  to  become  unstable,  the  even  mode  takes  over  at  roughly  X  =  82  cm  to  become 
the  most  unstable  mode.  The  growth  rate  of  the  secondary  instability  is  an  order  of  magnitude 
higher  than  that  of  the  nonlinear  Gortler  vortex. 

The  frequency  and  wavenumber  of  the  secondary  instability  corresponding  to  the  maximum 
growth  rate  at  various  streamwise  locations  are  given  in  Table  1.  We  see  that  the  frequencies  vary 
greatly  from  one  streamwise  station  to  the  next.  The  experimentally  observed  frequency  of  130  Hz 
may  not  be  that  of  the  most  unstable  wave  at  large  downstream  distances.  The  relatively  lower 
frequency  waves  (around  100  Hz)  become  unstable  first.  When  the  region  of  higher  frequency  waves 
is  reached,  the  laminar  basic  flow  state  may  have  already  been  destroyed  by  the  nonlinear  growth  of 
secondary  instabilities  and  the  higher  frequency  waves  may  not  have  a  chance  to  manifest 
themselves.  The  temporal  direct  numerical  simulation  of  Liu  &  Domaradzki  (1993)  has  a 
computational  box  which  restricts  the  maximuin  wavelength  to  2.0  cm.  According  to  Table  1,  the 
frequency  of  most  amplified  disturbance  is  just  above  160  Hz  (their  calculation  gives  about  200  Hz). 


Table  L  Frequency  and  wavelength  of  the  secondary  instability  modes 


Odd  Mode 

Even  Mode 

X(cm) 

/(Hz) 

Wavelength 

(cm) 

X(cm) 

/(Hz) 

Wavelength 

(cm) 

65.0 

67.1 

4.59 

79.8 

181.7 

1.66 

70.0 

87.8 

3.57 

82.3 

194.8 

1.56 

74.9 

108.8 

2.93 

84.9 

210.7 

1.47 

82.3 

141.2 

2.38 

89.7 

243.2 

1.34 

84.9 

161.6 

2.07 

94.6 

271.3 

1.27 

89.7 

202.1 

1.70 

94.6 

243.0 

1.47 

The  calculations  of  Hall  &  Horseman  (1991)  were  almost  exclusively  for  the  basic  flow  state  atX  = 
100  cm.  They  found  that  the  fastest  growing  wave  (odd  mode)  has  a  frequency  of  about  110  Hz  and  a 
wavelength  of  about  3  cm.  This  would  compare  well  with  our  results  at  about  X  =  76  cm.  This  can  be 
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explained  by  noting  the  fact  that  their  disturbed  flow  state  amplitude  is  apparently  lower  than  ours 
at  corresponding  streamwise  locations;  at  X  =  100  cm,  our  calculation  in  Figure  1  shows  a  thin  neck 
while  their  u(y,z)  plot  does  not  show  this  and  resembles  more  with  the  structure  shown  in  our 
Figure  2(b). 

We  now  consider  the  eigenfunctions  of  these  secondary  instability  modes  at  X  =  95  cm.  The 
wavelength  chosen  here  is  1.53  cm,  close  to  the  fastest  growing  odd  and  even  modes  at  this 
downstream  location.  In  addition  to  the  most  unstable  odd  and  even  modes,  a  second  even  mode  is 
also  considered.  The  contours  of  velocity  eigenfunction  |i2g|  atX  =  95  cm  are  shown  as  solid  lines  in 
Figure  10,  together  with  basic  flow  state  u(y,z)  as  dashed  lines  in  the  background.  The 
eigenfunctions  are  normalized  so  that  the  maximum  |ug|  has  an  amplitude  of  unity.  The  contours 
plotted  are  from  0.1  to  0.9  in  intervals  of  0.1,  and  u{y,z)  contours  plotted  are  from  0.1  to  0.9  in 
intervals  of  0.1  U^.  One  feature  to  notice  is  that  the  phase  speeds,  c^,  of  the  modes  shown  in  Figure 
10  are  all  close  to  0.7  U^,  and  the  amplitudes  of  these  modes  are  concentrated  in  the  neighborhood  of 
the  manifold  u{y,z)  =  0.7  U^.  This  manifold  would  be  the  critical  layer  in  the  case  of  neutral 
stability. 

The  contours  shown  in  Figure  10  bear  some  similarity  to  those  obtained  by  Hall  &  Horseman 
(1991)  despite  the  difference  between  the  basic  flow  states  in  their  work  and  the  present  work.  The 
odd  mode  has  two  dominant  peaks,  one  on  either  side  of  the  peak  plane.  In  the  case  of  the  even 
modes,  the  second  even  mode  has  three  dominant  peaks  similar  to  that  shown  in  Figure  6  (a)  of  Hall 
&  Horseman  (1991)  for  the  only  even  mode  they  analyzed  in  their  work.  The  most  unstable  single¬ 
peak  even  mode  was  not  mentioned  in  Hall  &  Horseman  (1991).  Considering  the  fact  that  the 
growth  rate  of  the  first  odd  mode  is  about  1.7  times  that  of  the  second  even  mode  and  that,  in  the 
work  of  Hall  &  Horseman  (1991),  the  odd  mode  grows  almost  twice  as  fast  as  the  even  mode,  we  are 
led  to  believe  that  the  three-peak  even  mode  analyzed  by  Hall  &  Horseman  (1991)  was  actually  the 
second  unstable  even  mode,  and  the  first  unstable  even  mode  was  missed  because  they  used  a 
shooting  technique  to  compute  the  eigenvalues.  In  our  study,  we  use  a  global  method  which  finds  all 
the  eigenvalues  of  the  discretized  problem. 

The  contours  of  (u,(  for  the  odd  mode  bears  striking  similarity  to  the  streamwise  rms 
fluctuations  shown  in  Figure  16  of  the  temporal  simulation  of  Liu  &  Domaradzki  (1993),  suggesting 
that  this  mode  indeed  plays  an  important  role  in  the  break-up  of  Gdrtler  vortices. 

From  the  eigenfunctions,  we  see  a  definite  relationship  between  instability  and  inflection  in  the 
velocity  profiles.  Figure  11  shows  contours  of  vertical  and  horizontal  velocity  gradients  Uy  and  u^. 
The  inflection  in  the  velocity  profiles  occurs  at  points  where  the  velocity  gradient  is  maximum.  The 
amplitudes  of  the  eigenfunctions  concentrate  near  the  regions  of  maximum  velocity  gradient. 
Furthermore,  the  most  unstable  even  mode  is  clearly  associated  with  the  vertical  velocity  gradient. 
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the  most  unstable  odd  mode  is  associated  with  the  horizontal  velocity  gradient  and  the  second  most 
unstable  even  mode  is  associated  with  both  gradients. 

We  have  analyzed  eigenmodes  which  are  unstable  at  different  streamwise  locations  if  the  base 
flow  of  Gdrtler  vortex  is  free  of  unsteady  perturbations  of  large  amplitude.  It  will  be  shown  in  the 
next  section  that  nonlinearity  will  cause  the  rms  amplitude  of  the  secondary  instability  to  be 
modified. 

Finally,  we  point  out  that  the  relative  importance  of  odd  and  even  modes  may  depend  upon  the 
flow  parameters  (such  as  Gdrtler  number,  wave  number,  etc.)  which  govern  the  evolution  of  the 
steady  vortices.  However,  the  general  "shape"  of  the  eigenfunctions  and  their  association  with  the 
vertical  or  spanwise  velocity  gradients  should  remain  the  same. 

5.  Nonlinear  Development  of  Unsteady  Disturbances 

We  now  solve  Eqs.  (2.7-2.8)  to  study  nonlinear  evolution  of  the  steady  as  well  as  unsteady 
disturbances.  The  spatial  secondary  instability  computations  for  the  even  and  odd  modes  are  started 
at  streamwise  locations  where  the  respective  modes  are  moderately  unstable.  Since  exact  initial 
conditions  are  diflicult  to  obtain,  we  use  the  eigenfunctions  obtained  from  inviscid  linear  secondary 
instability  analysis  to  approximate  these  conditions.  Calculations  show  that  transients  decay  very 
fast.  The  initial  amplitude  assigned  to  these  disturbances  is  small  enough  (of  in  u,  for  both 

odd  and  even  modes)  to  ensure  that  the  initial  evolution  of  the  secondary  instability  is  linear. 

We  first  discuss  the  computations  for  the  odd  mode.  These  calculations  are  performed  in  the 
following  way.  We  start  the  calculations  at  X  =  10cm  for  only  the  steady  disturbances.  These 
calculations  are  carried  up  to  X  =  75  cm  where  the  unsteady  odd  mode  disturbances  are  introduced. 
The  frequency  of  the  fundamental  secondary  instability  is  chosen  to  be  110  Hz,  very  close  to  the  most 
unstable  odd  mode  at  that  location  (109  Hz).  The  number  of  Fourier  modes  is  11  in  the  spanwise 
direction  (-10  re  S  10)  and  8  in  time  (-7  <  /re  <  7).  Figure  12  shows  contours  of  streamwise  rms 
fluctuations  at  four  streamwise  locations  downstream  of  the  starting  location.  Initially,  the  shape  of 
amplitude  distribution  very  much  resembles  the  local  eigenfunctions  analyzed  in  the  last  section. 
Later  on,  at  larger  x,  nonlinearity  causes  the  amplitude  distribution  to  become  fatter,  beginning  to 
fill  up  the  y-z  plane.  The  maximum  amplitude  reaches  about  20  percent  of  the  free  stream  velocity, 
.  There  are  two  regions  of  high  amplitude:  one  near  the  wall  and  the  other  away  from  the  wall. 
Initially  the  region  near  the  boundary-layer  edge  has  higher  amplitude,  but  as  the  disturbances 
evolve  downstream  the  near-wall  region  attains  higher  amplitude.  The  contours  of  streamwise  mean 
flow  (time  averaged  flow)  at  corresponding  streamwise  locations  are  shown  in  Figure  13.  The 
changes  in  the  "mushroom"  structure  from  that  without  the  high  frequency  oscillations  far 
downstream  are  profound  (compare  with  Figure  1(c)).  The  top  of  the  "mushroom"  becomes  flatter 
due  to  the  fact  that  strong  unsteady  oscillations  smear  out  the  differences  in  velocity  gradients.  As 
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the  oscillations  get  stronger,  the  "thin  neck"  region  of  the  "mushroom",  where  low  speed  fluid  lifts  up, 
becomes  conflned  to  the  neighborhood  close  to  the  wall.  The  experiments  of  SB  show  similar 
structure  of  the  vortices  prior  to  breakdown,  except  that  in  their  experiment  vortices  tend  tu  loose 
S3mimetry  about  the  peak  plane.  Some  of  our  calculations  suggest  that  such  loss  of  S3rmmetry  is 
linked  to  the  unsymmetric  initial  disturbance  Held. 

We  now  carry  out  computations  for  the  even  mode,  which  is  introduced  at  the  streamwise 
location  X  =  82  cm.  The  frequency  of  the  fundamental  secondary  instability  is  taken  to  be  195  Hz, 
which  is  the  frequency  of  the  most  unstable  even  mode  at  that  location.  The  number  of  Fourier 
modes  is  the  same  as  in  the  odd  mode  case.  In  the  early  linear  stage,  the  distribution  of  amplitude 
of  the  streamwise  rms  fluctuation  is  concentrated  in  the  region  away  from  the  wall.  Nonlinear 
interaction  bring  the  perturbation  level  near  the  wall.  The  rms  values  of  streamwise  fluctuations 
and  the  mean  flow  are  shown  in  Figures  14  and  Figure  15,  respectively. 

In  their  numerical  simulation,  Liu  &  Domaradzki  (1993)  observed  that  the  oscillation  frequency 
of  the  vertical  velocity  is  twice  that  of  the  spanwise  velocity  in  the  low  speed  region,  while  the  two 
frequencies  are  the  same  away  from  that  region.  This  is  also  found  to  be  the  case  in  the  present 
computations  for  the  odd  mode.  Figure  16  shows  the  vertical  and  spanwise  velocity  fluctuations  at 
two  fixed  locations  in  space,  one  in  the  low  speed  region,  the  other  away  from  it.  Counting  the 
number  of  dominant  peaks,  we  can  observe  the  difference  in  frequency  in  the  low  speed  region.  This 
phenomenon  can  be  explained  by  considering  the  spatial  symmetry  preservation  property  of  the 
fundamental  secondary  instability  and  its  harmonics.  Initially,  the  only  mode  present  is  the 
fundamental  with  frequency,  say,  f.  In  the  odd  mode  case,  the  vertical  velocity  fluctuation  for  the 
fundamental  mode  is  odd  with  respect  to  the  line  of  symmetry  of  the  background  "mushroom" 
structure,  while  the  spanwise  velocity  fluctuation  is  even.  It  can  be  mathematically  verified  that  the 
symmetry  properties  of  the  harmonics  obey  the  following  rule:  for  the  vertical  velocity,  modes  with 
frequencies  2f,  4f,  6f,  etc.  are  even,  while  modes  with  frequencies  3f,  5f,  If,  etc.  are  odd;  for  the 
spanwise  velocity,  modes  with  frequencies  2f,  4f,  6f,  etc.  are  odd,  while  modes  with  frequencies  3f,  5f, 
If,  etc.  are  even.  These  symmetry  properties  are  preserved  as  the  flow  develops  downstream. 
Suppose  we  place  a  velocity  probe  somewhere  along  the  line  of  symmetry  of  the  background 
"mushroom"  structure,  we  cannot  detect  the  amplitude  of  the  fundamental  mode  of  the  vertical 
velocity  since  it  is  odd  and,  therefore,  has  zero-amplitude  there.  The  lowest  frequency  mode  that  can 
be  detected  is  the  2f  mode  of  the  vertical  velocity.  The  lowest  frequency  mode  of  the  spanwise 
velocity  that  can  be  detected  is  the  fundamental  mode.  Hence,  it  appears  that,  in  the  peak  region, 
the  frequency  of  the  vertical  velocity  is  twice  that  of  the  spanwise  velocity.  Once  we  move  away  from 
the  peak  region,  the  fundamental  modes  of  both  the  vertical  and  the  spanwise  velocities  can  be 
detected,  therefore,  the  same  frequencies  are  observed. 
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We  define  the  energy  associated  with  frequency,  mf,  as 


=  X  +K«f 


As-mO 

for  m  >  0,  and 


A=1 


j(l«0«f  +|«0nr  +Knfy>'  +|j(l“0of +I‘"00r)cf3' 


(5.1a) 


(5.1b) 


and  plot  it  in  Fig.  17  as  a  function  of  down  stream  distance  for  both  the  sinuous  and  the  varicose 
secondary  instability.  The  change  in  the  zero-frequency  mode  is  relatively  slow.  The  energy  curves 
of  the  fundamental-frequency  modes  for  both  type  of  secondary  instabilities  are  roughly  straight 
lines  over  a  considerable  range  of  streamwise  distance  before  the  energies  in  the  harmonics  become 
significant,  indicating  the  rapid  exponential  growth. 

We  plot  the  instantaneous  streamwise  velocity  contours  in  the  x-z  plane  at  y  =  1.08  cm  (Figure 
18).  We  see  that  the  odd  mode  perturbs  the  Gdrtler  vortices  in  a  wavy  (or  sinuous  )  manner,  while 
the  even  mode  breaks  up  the  otherwise  straight  contours  of  Gdrtler  vortices  into  series  of  knotty 
structures  associated  with  the  "horse-shoe"  vortex  mode  of  breakdown.  At  large  downstream 
distances,  more  and  more  small  scale  structures  begin  to  appear,  as  the  flow  heads  for  transition  to 
turbulence. 

In  Fig.  19,  we  plot  the  instantaneous  streamwise  velocity  (snapshots  at  three  fixed  times)  in  the 
y-z  plane  at  x  =  101  cm  for  the  sinuous  mode.  Comparison  with  Fig.  1  clearly  shows  the  effect  of 
unsteady  disturbances  on  the  Gdrtler  vortex.  The  unsteady  motion  computed  here  is  qualitatively 
similar  to  that  observed  in  the  experiments  (e.g.,  Peerhossaini  &  Wesfreid  (1988). 

Finally,  we  point  out  that,  as  in  the  case  of  the  steady  Gdrtler  vortex  discussed  earlier, 
quantitative  one-to-one  comparison  with  the  experiment  of  SB  is  made  difficult  by  the  uncertainty  in 
the  initial  amplitude  of  the  secondary  instability  and  the  spanwise  variation  of  (jdrtler  vortex 
wavelength  in  the  experiment.  For  example,  the  computed  isocontours  of  streamwise  velocity  for  the 
odd  mode  at  x  =  102.3cm  shown  in  Fig.  13  agree  much  better  with  results  of  SB  at  110  cm  than  at 
100  cm. 


6.  Conclusions 

Computation  of  the  nonlinear  development  of  steady  Gdrtler  vortices,  their  stability 
characteristics  with  respect  to  high  frequency  secondary  disturbances  and  the  nonlinear  spatial 
development  of  the  secondary  instabilities  are  carried  out.  Qualitative  agreement  with  the 
experiment  of  Swearingen  &  Blackwelder  (1987)  (SB)  is  obtained.  It  is  found  that  the  wall  shear  in 
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the  peak  plane  begins  to  increase  at  a  sufllciently  large  downstream  distance  due  to  the  nonlinear 
interaction  of  the  stationary  modes  even  before  the  unsteady  oscillations  become  strong.  The  cause 
of  this  phenomenon  is  the  fact  that  the  mean  flow  correction  mode  ultimately  becomes  the  dominant 
mode  and  overcomes  the  effort  of  the  Gdrtler  vortices  to  slow  down  the  flow  in  the  peak  region.  The 
energy  in  each  stationary  mode  eventually  approaches  saturation.  The  computed  "mushroom" 
structures  bear  strong  resemblance  to  those  obtained  in  the  experiment  of  SB. 

For  this  particular  basic  flow  state  which  is  set  up  due  to  the  Gdrtler  vortices,  the  temporal 
secondary  instability  analysis  is  carried  out  using  a  2-D  eigenvalue  approach  associated  with  the 
governing  partial  differential  equations.  For  the  conditions  of  the  experiment  of  SB,  the  odd  mode  of 
secondary  instability  begins  to  show  up  at  approximately  X  =  60  cm,  the  even  mode  becomes 
unstable  later  at  approximately  70  cm  <  X  <75  cm.  At  about  X=82  cm,  the  even  mode  becomes  more 
unstable  than  the  odd  mode.  Comparisons  of  the  amplitude  distributions  of  eigenfunctions  with  the 
distributions  of  vertical  and  spanwise  shear  gradients  clearly  indicate  the  close  association  of  the 
most  unstable  even  mode  with  the  vertical  shear  and  the  most  unstable  odd  mode  with  the  spanwise 
shear. 

The  nonlinear,  spatial  development  of  the  odd  and  even  modes  of  secondary  instability  is 
computed  using  the  PSE  method.  The  odd  and  even  modes  give  rise  to  the  sinuous  instability  and 
the  varicose  instability,  respectively.  Either  mode  can  lead  to  the  breakdown  of  Gbrtler  vortices. 
The  two-fold  difference  between  the  frequency  of  the  vertical  velocity  and  that  of  the  spanwise 
velocity  oscillations  found  by  Liu  &  Domaradzki  (1993)  is  simply  due  to  the  fact  that  the 
fundamental  mode  of  the  vertical  velocity  has  zero  amplitude  at  the  line  of  symmetry  of  the 
"mushroom"  for  the  sinuous  instability  and  the  high  frequency  in  the  vertical  velocity  they  detected 
was,  in  fact,  for  the  harmonic.  Prior  to  breakdown  to  turbulence,  the  nonlinear  interaction  among 
the  steady  and  unsteady  modes  eventually  make  Gdrtler  vortices  to  oscillate  sinuously  in  the  plane 
parallel  to  the  plate  or  cause  them  to  develop  horse-shoe  type  structures  which  travel  downstream. 
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Fig.  4.  Evolution  of  disturbance  energy  for  various  spanwise  Fourier  modes. 
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Fig.  5.  Normal  wall  shear  at  the  peak  and  the  valley.  Symbols  represent  the  experimental  data  of  SB. 
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Fig.  7.  Evolution  of  displacement  thickness  at  the  peak  and  the  valley.  S3nnbols  represent  the  experimental  data  of  SB. 
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Fig.  9.  Maximum  spatial  growth  rate  of  the  odd  and  even  modes  at  various  streamwise  locations,  compared  with  the  growth  rate 
Giirtler  vortex. 
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Variation  of  time-averaged  streamwise  velocity  in  the  Y-Z  plane  atX  =  102.3  cm.  Odd  mode. 
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Variation  of  rms  streamwise  velocity  fluctuations  in  Y-Z  plane  for  the  even  mode  (a)X=  87.8  cm;  (b)X=  100.4  cm;  (c)X  =  102.3  cm; 
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Fig.  16.  Vertical  and  spanwise  velocity  fluctuations  with  time  at  X  =  105  cm  and  Y  =  1.29  cm.  Upper  graph  for  fluctuations  at  the  peak. 
Lower  graph  for  fluctuations  at  Z  =  XJA  away  from  the  peak. 
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Fig.  18.  Instantaneous  streamwise  velocity  in  the  X-Z  plane  at  F  =  1.08  cm.  Upper  half  —  sinuous  mode  (odd),  Lower  half 
mode  (even). 
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Fig.  19.  Instantaneous  streamwise  velocity  in  they-z  plane  for  the  sinuous  mode  (/=  110  Hz):  top, 
t  =  0;  middle,  t  =  ^rdZar,  bottom,  t  =  4«/3g}. 
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PARTB 


Crossflow  Disturbances  in  Three-Dimensional  Boundary 
Layers:  Nonlinear  Development,  Wave  Interaction  and 

Secondary  Instability 
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Abstract 

Nonlinear  stability  of  a  model  swept-wing  boundary  layer,  subject  to  crossflow  instability,  is 
investigated  by  numerically  solving  the  governing  partial  differential  equations.  The  three- 
dimensional  boundary  layer  is  unstable  to  both  stationary  and  traveling  crossflow  disturbances. 
Nonlinear  calculations  have  been  carried  out  for  stationary  vortices  and  the  computed  wall  vorticity 
pattern  results  in  streamwise  streaks  which  resemble  quite  well  with  the  surface  oil-flow 
visualizations  in  swept-wing  experiments.  Other  features  of  the  stationary  vortex  development  (half¬ 
mushroom  structure,  inflected  velocity  profiles,  vortex  doubUng,  etc.)  are  also  captured  in  these 
calculations.  Nonlinear  interaction  of  the  stationary  and  traveling  waves  is  also  studied.  When 
initial  amplitude  of  the  stationary  vortex  is  large  as  compared  to  the  traveling  mode,  the  stationary 
vortex  dominates  most  of  the  downstream  development.  When  the  two  modes  have  the  same  initial 
amplitude,  the  traveling  mode  dominates  the  downstream  development  owing  to  its  higher  growth 
rate.  It  is  also  found  that,  prior  to  laminar/turbulent  transition,  the  three-dimensional  boundary 
layer  is  subject  to  a  high  frequen<y  secondary  instability  which  is  in  agreement  with  the  experiments 
of  Poll  (1985)  and  Kohama,  Saric  &  Hoos  (1991).  The  fi*equency  of  this  secondary  instability,  which 
resides  on  top  of  the  stationary  crossflow  vortex,  is  an  order  of  magnitude  higher  than  the  frequency 
of  the  most  amplified  traveling  crossflow  mode. 


1.  Introduction 

In  swept-wing  flows,  chordwise  pressure-gradient  near  the  leading  edge  causes  inviscid  stream¬ 
lines  to  be  curved  in  the  planes  parallel  to  the  wing  surface.  Associated  with  this  streamline 
curvature  is  a  pressure  gradient  which  acts  in  a  direction  normal  to  the  streamlines  and  introduces  a 
secondary  flow  within  the  boundary-layer.  This  secondsuy  flow,  commonly  known  as  crossflow,  is 
subject  to  inviscid  instability  due  to  the  presence  of  an  inflection  point  (Gregory,  Stuart  &  Walker 
1955)  and  is  the  main  cause  of  transition  in  swept-wing  flows.  Thus,  this  problem  is  not  only  of 
fundamental  importance  in  fluid  mechanics  but  also  of  prime  significance  in  laminar  flow  control 
(LFC)  design  of  swept-wings. 

Crossflow  instability  often  results  in  the  formation  of  stationary  corotating  vortices  commonly 
called  crossflow  vortices.  This  phenomenon  is  observed  in  swept-wing  boundary  layers  as  well  as  in 
other  geometries  such  as  rotating  disks  and  cones.  How  the  stationary  crossflow  vortices  lead  to 
turbulence  remains  unanswered.  Traveling  crossflow  disturbances  are  also  possible  and  the  role  of 
traveling  vs.  stationaiy  disturbances  is  a  question  which  needs  to  be  investigated.  Another  problem 
which  is  of  interest  in  swept-wing  flows  is  the  possibility  of  interaction  between  the  inviscid  crossflow 
disturbances  and  viscous  streamwise  instability.  Crossflow  disturbances  are  amplified  in  the 
negative  pressure  rise  region  near  the  wing  leading  edge,  while  Tollmien-Schlichting  (TS)  waves 
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(viscous  instability  of  streamwise  profiles)  are  amplified  in  the  flat  pressure  region  of  the  wing 
midchord.  The  possible  interaction  of  these  two  types  of  disturbances  may  be  quite  significant  in  the 
successful  design  of  LFC  wings. 

Experimental  investigations  into  the  nature  of  the  swept-wing  boundary-layer  instability,  at  its 
linear  and  nonlinear  stage,  have  been  carried  out  by  Bippes  and  coworkers  (see,  Bippes  (1991), 
Miiller  &  Bippes  (1988),  Muller  (1989))  at  DLR,  by  Saric  and  coworkers  (see,  Dagenhsirt  et  al.  (1989), 
Saric,  Oaganhart  &  Mousseux  (1989),  Kohama,  Saric  &  Hoos  (1991))  at  Arizona  State  University  and 
by  Amal  and  coworkers  (see  Amal  &  Juillen  (1987))  at  ONERA/CERT.  Experiments  at  DLR  were 
performed  on  a  swept-plate  model  with  an  imposed  pressure  gradient.  A  displacement  body  above 
the  plate  was  used  to  generate  the  Cp  distribution  which  varied  almost  linearly  with  chordwise 
distance.  Saric  and  Amal  used  infinite-swept  aerofoils  in  their  low-speed  experiments. 

Both  stationary  and  traveling  crossflow  disturbances  were  observed  in  these  experiments,  as 
well  as  in  the  experiment  of  Poll  (1985)  on  a  swept  cylinder.  Muller  &  Bippes  (1988)  found  that  the 
stationary  vortices,  as  well  as  traveling  disturbances,  reached  nonlinear  saturation  in  their 
experiment  However,  they  did  not  notice  any  explosive  secondary  instability  leading  to  transition. 
On  the  other  hand,  Kohama,  Saric  &  Hoos  (1991)  observed  a  high  frequency  secondary  instability 
prior  to  transition  in  their  swept-wing  experiment  where  pressure  gradient  remained  favorable 
ruling  out  any  possibility  of  TS  wave  amplification.  The  frequency  of  this  secondary  instability  was 
an  order  of  magnitude  higher  than  the  frequency  of  the  most  amplified  traveling  disturbance  given 
by  the  linear  theory.  They  concluded  that,  even  though  the  traveling  crossflow  disturbances  are 
observed,  the  transition  process  in  this  three-dimensional  boundary  layer  is  dominated  by  the 
stationary  vortices  and  the  associated  secondary  instability.  Poll  (1985)  had  also  observed  a  high 
frequency  disturbance  in  his  swept-cylinder  experiment. 

Miiller  &  Bippes  (1988)  also  studied  the  effect  of  free-stream  turbulence  on  the  instability 
behavior  in  their  experiment.  They  found  that  at  "low"  levels  (.05  percent)  of  free  stream  turbulence, 
stationary  disturbances  amplified  to  large  amplitudes  but  these  large  amplitudes  of  the  stationary 
vortices  did  not  necessarily  lead  to  early  transition.  The  experiments  performed  in  wind  tunnels 
with  higher  turbulence  levels  (.15  and  .3  percent)  showed  weaker  growth  of  stationary  disturbances 
but  earlier  transition  due  to  stronger  traveling  disturbances.  They  concluded  that  traveling  waves, 
and  not  the  stationary  vortices,  play  the  major  role  in  the  transition  process.  Their  experimental 
results  also  seem  to  suggest  an  early  nonlinear  interaction  between  stationary  and  traveling 
crossflow  disturbances. 

Theoretical  investigations  into  linear  and  nonlinear  stability  of  three-dimensional  boundary  lay¬ 
ers  have  been  carried  out  by  Balachandar,  Streett  &  Malik  (1992),  Fischer  &  Dallmann  (1991),  Malik 
(1986),  Meyer  &  Kleiser  (1988)  and  Reed  (1987).  Fischer  and  Dallman  used  secondary  instability 
theory  and  Meyer  and  Kleiser  used  direct  simulation  of  Navier-Stokes  equations  to  study  the  swept- 
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plate  experiment  of  Muller  &  Bippes  (1988).  Fischer  and  Dallmann  argued  that  the  traveling 
disturbances  observed  in  the  DLR  experiment  are  secondary  disturbances  of  the  mean  flow  modu¬ 
lated  by  the  stationary  vortices  and  should  not  be  thought  of  as  the  primary  instability  of  the  three- 
dimensional  boimdary-layer  flow.  Direct  numerical  simulation  of  Meyer  and  Kleiser  found  nonlinear 
equilibrium  states  for  stationary  as  well  as  traveling  disturbances,  in  agreement  with  the  DLR 
experiment.  Similar  equilibrium  states  for  stationary  vortices  were  computed  by  Malik  (1986)  in 
rotating-disk  boundary  layer.  Balachandar,  Street!  &  Malik  (1992)  performed  a  secondary 
instability  analysis  of  the  rotating-disk  boundary  layer  where  the  stationary  vortices  constituted  the 
primary  instability.  They  were  able  to  find  a  high-frequency  secondary  instability  similar  to  the  one 
observed  by  Kohama  (1984,1987)  in  a  rotating-disk  boundary  layer. 

All  these  experimental  and  theoretical  investigations  consider  a  class  of  mean  flows  which  is 
only  unstable  to  inflectional  crossflow  disturbances  and  do  not  support  TS  wave  amplification. 
Results  from  various  experiments  appear  to  suf^est  that  for  this  class  of  flow,  there  are  at  least  two 
possible  scenarios  for  transition.  If  free-stream  turbulence  level  is  very  small,  i.e.,  the  initial 
amplitude  of  the  nonstationary  disturbances  is  small  relative  to  the  stationary  disturbances,  which 
most  certainly  are  introduced  at  local  surface  imperfections,  then  stationary  disturbances  dominate 
the  initial  stage  of  the  disturbance  growth  leading  to  a  high-frequen<y  secondary  instability  resulting 
in  final  breakdown.  When  the  initial  amplitude  of  the  traveling  modes  is  not  small,  nonlinear 
interaction  between  these  traveling  modes  and  stationary  vortices  is  present  and  the  character  of  the 
final  breakdown  is  influenced  by  the  relative  amplitudes  of  the  stationary  vortices  and  the  traveling 
modes.  The  other  class  of  flow  where  TS  waves  could  amplify  is  also  of  technological  importance  but 
has  not  been  studied  either  experimentally  or  theoretically. 

The  objective  of  this  research  is  to  study  various  wave-interaction  mechanisms  and  laminar-flow 
breakdown  in  three-dimensional  boundaiy  layers.  Previous  linear  and  nonlinear  theoretical  investi¬ 
gations  have  been  performed  by  using  parallel-flow  approximation  and  have  been  local  in  nature. 
This  study  includes  nonparallel  effects  and  sets  up  the  problem  within  the  framework  of  nonlinear 
parabolized  stability  equations  (PSE).  Intermodal  interaction  and  the  effect  of  initial  conditions  can 
also  be  studied  by  using  this  approach.  Basic  insight  into  the  physical  mechanisms  involved  in 
swept-wing  flow  transition  can  be  achieved  by  considering  simple  model  flows.  One  such  flow  is  the 
swept  Hiemenz  flow  in  which  the  interaction  of  stationary  and  traveling  crossflow  disturbances  can 
be  studied.  In  this  paper  we  study  linear  and  nonlinear  crossflow  disturbances  as  well  as  the 
interaction  between  stationary  and  traveling  modes.  We  also  study  secondary  instability  of  the 
three-dimensional  mean  flow  modulated  by  the  stationary  vortices.  Section  2  describes  the  basic 
flow  for  the  swept  Hiemenz  problem  and  the  associated  PSE  analysis  is  given  in  §  3.  The  results  for 
linear  and  nonlinear  stabiHty  analysis  and  wave  interactions  are  given  in  §  4.  Section  5  describes  the 
results  from  secondary  instability  analysis  and  the  conclusions  are  given  in  §  6. 
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2.  The  Swept  Hiemenz  Problem 

The  flow  past  a  circular  cylinder,  outside  the  viscous  boundary-layer,  can  be  represented  as 

=CX*  +CiX*^  +02***  +  .  .  .  (2.1) 

where  is  the  velocity  along  the  coordinate  **  and  c,Ci,  etc.  are  constants.  In  the  two-dimensional 
stagnation-point  flow,  only  the  first  term  in  the  series  (2.1)  is  retained  and,  hence,  the  velocity 
increases  linearly  with  distance  x*  ,  i.e., 

f/.  =cx*  (2.2) 

If  we  consider  the  Cartesian  coordinate  system  x*,  y*,  z*,  then  (2.2)  gives  the  far-field  (y* 
solution  of  the  impinging  flow  on  a  plate  along  x*  which  we  define  here  by  y*  =  0.  The  associated 
viscous  problem  was  first  investigated  by  Hiemenz  who  found  an  exact  solution  which  is  named  after 
him.  This  stagnation  flow  was  found  to  be  stable  to  infinitesimally  small  disturbances  propagating 
along  z*  by  Wilson  &  Gladwell  (1978). 

The  swept  Hiemenz  problem  is  constructed  by  introducing  a  velocity  component  along  the  z* 
axis  which  amounts  to  changing  the  inclination  of  the  impinging  stream  with  respect  to  z*.  The  flow 
is  symmetric  about  the  line  x*  =  0  which  is  called  the  attachment-line.  Linear  and  nonlinear 
stability  of  the  attachment- line  boundaiy-layer  has  been  studied  by  Hall,  Malik  &  Poll  (1984)  and 
Spalart  (1988).  In  this  paper  we  study  the  stability  of  this  flow  for  x*  >  0  as  was  recently  done  by 
Spalart  (1989)  using  full  Navier-Stokes  equations. 

2.1.  The  basic  flow 

We  consider  the  flow  of  a  viscous  incompressible  fluid  of  kinematic  viscosity  v.  Let  t  =  Vv/c  be 
a  typical  thickness  of  the  boundary  layer  which  is  used  here  as  the  length  scale.  We  note  that  t  is 
independent  of  x*.  Thus,  we  have  the  scaled  coordinates  x,y,z  given  as 

♦  «  • 

/  \  y  ^  \ 

(x,y,z)-(  ^  ^  • 


We  also  define  two  Reynolds  numbers  R  and  R  where 

II 

(2.3) 

V 

(2.4) 

From  (2.2)  and  (2.3),  it  follows  that 

R  =  —  =  x 

t 

(2.5) 

The  local  angle  of  the  inviscid  streamline  6,  with  respect  to  the  x-axis,  is  given  as 
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w 

0  =  atani  1  =  atan 


(2.6) 


We  now  look  for  a  solution  to  the  Navier-Stokes  equations  which  satisfies  the  following 
conditions 

u'  =v*  =w*  =0,  y*  =0  (2.7) 


(2.8) 


where  u*,  v*,  w*  are  velocity  components  in  the  x*,  y*,  z*  directions,  respectively.  It  is  convenient  to 


define  a  stream  function  0  so  that 


•  d0  •  d0 


and 


<P  =  X*  4^ f{y) 

If  we  use  W„  as  the  velocity  scale,  then 


Similarly, 


-  “  ^  f’{  \ 


—  tv  ,  . 


where  f  andg  are  governed  by  the  ordinary  differential  equations 

/•"'  +  /'/-"  +  (l-r^)  =  0 

g"  +  fg'  =  0 

where  primes  denote  differentiation  with  respect  toy. 

The  mean  flow  derivatives  needed  in  the  stability  analysis  below  can  be  written  as 

u  =  —  f"  u  =-ir" 

iVy  =g',  Wy  =g" 

-  1 

“x  =■=/ 

*  R 


(2.9) 

(2.10) 

(2.11) 

(2.12) 

(2.13) 
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Numerical  solution  of  (2.12)  and  (2.13)  thus  yields  the  mean  flow  and  its  derivatives  without  any 
additional  approximation  and  the  boundary-layer  thickness  does  not  vary  with  x. 

3.  PSE  Analysis  for  3D  Boundary  Layers 

Parabolized  stability  equations  (PSE)  for  linear  and  nonlinear  disturbances  in  two-dimensional 
boundary  layers  have  been  used  by  Herbert  (1991)  and  Bertolotti,  Herbert  &  Spalart  (1992)  for 
incompressible  flow  where  they  used  streamfunction  formulation  of  the  governing  equations.  In  the 
present  three-dimensional  (3D)  boundary-layer  study,  we  follow  the  work  of  Chang  et  al.  (1991)  for 
compressible  flow  and  formulate  the  incompressible  stability  problem  using  primitive  variables  in 
Cartesian  coordinates  x,y,  and  z.  The  basic  flow  is  perturbed  by  fluctuations  in  the  flow,  i.e.,  the 
total  field  can  be  decomposed  into  a  mean  value  (solution  of  (2.12-2.13))  and  a  perturbation  quantity 

u  =  u+u,  v  =  v+v,  w  =  w+w,  p  =  p  +  p  (3.1) 


where  p  is  the  pressure.  Substituting  (3.1)  into  the  incompressible  Navier-Stokes  equations  and 
subtracting  from  it  the  steady  mean  flow,  we  obtain  the  nonlinear  disturbance  equations  as 


dt  dx  dy  dz 


(3.2) 


where  the  left  hand  side  contains  only  linear  operators  operating  on  the  disturbance  vector 
0  =  (u,o,io,p)  and  the  right-hand-side  forcing  vector  F  is  due  to  non-linear  interaction  and  includes 
all  non-linear  terms  associated  with  the  disturbances.  The  right  hand  side  is  given  as 


dx  dy 


^  dz  ' 


(3.3) 


In  the  above,  r  is  the  diagonal  matrix  [1,1,1,0]  while  A,  B,  C  are  given  as 


A  = 


and  A,  B,  C  are  similar  to  A,  B,  C  except  that  quantities  with  over  bar  are  replaced  with  ~  and  all 
ones  are  dropped.  The  coefficient  matrices  D,  Ey,  are  given  as 

Uy  0  0 


u 

0 

0 

1 

V 

0 

0 

O' 

w 

0 

0 

O' 

0 

u 

0 

0 

0 

V 

0 

1 

,  c  = 

0 

w 

0 

0 

0 

0 

u 

0 

0 

0 

V 

0 

0 

0 

w 

1 

.1 

0 

0 

0, 

_0 

1 

0 

0_ 

0 

0 

1 

0 

D  = 


Oj  Oy  0  0 


Wy  0  0 


0  0  0  0 
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R 


E,=Ey  =  E,= 


0 


0 

0 


0 

J_ 

R 

0 

0 


0  0 
0  0 


0  Oj 


We  assume  that  the  given  disturbance  is  periodic  in  time  and  in  the  spanwise  direction;  thus, 
the  disturbance  function  ^  can  be  expressed  by  the  following  Fourier  series 


0=  X  X 


(3.4) 


Here,  the  frequency  /a  and  wave  number  P  are  chosen  such  that  the  longest  period  and  wave  length 
are  2t^a)  and  2n/p  in  the  temporal  and  spanwise  domains,  respectively.  For  most  stability  problems 
of  interest,  it  is  sufficient  to  truncate  (3.4)  to  only  a  finite  number  of  modes 

u-i  N-i  (3.5) 

0=  X  X 

m=-Afn=-JV 


where  Af  and  N  represent  one-half  the  number  of  modes  kept  in  the  truncated  Fourier  series. 
Substituting  (3.5)  in  (3.2)  we  obtain  governing  equations  for  Xmn  which  are  elliptic.  In  order  to 
facilitate  the  solution  of  these  equations  we  decompose  the  disturbance  into  a  fast  varying  wave-like 
part  and  a  slowly  varying  shape  function  and  write  Xna 

Xmn(x,y)='i'nn(x,y)iirnn(x)  (3.6a) 


i  f 

Ji^n(x)  =  e  >0 


(3.6b) 


where  is  the  shape  function  etc.)  for  the  Fourier  mode  (monnp)  and  a„„  is  the 

associated  streamwise  (complex)  wave  number.  With  a  proper  choice  of  a„„  in  (3.6b),  the 
arbitrariness  in  (3.6a)  can  be  removed  and  the  equations  for  can  be  parabolized.  In  other  words, 
a^n  is  chosen  such  that  variation  of  with  x  is  minimized  which  allows  the  approximation 
-  0-  Th®  parabolized  stability  equations  (PSE)  for  the  shape  function  of  a  single  Fourier 
mode  (m,n)  can  be  written  as 

G.nn'l'r.n  ^  ^  =  E^  (3.7) 

where  matrices  E„„,  and  B„,„  are  given  by 

G„„  =  -imoir  +  ia„„A  +  inpC  +  D-Ej^i^^~  aL  j  + 
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A„„=A-2ia„„E, 


B„.n=B 


The  non-linear  forcing  function  is  the  Fourier  component  of  the  total  forcing,  F,  and  can  be 
evaluated  by  the  Fourier  series  expansion 

W-l  iV-l 

F(x.y,z,t)=  (3.8) 

m--M n=-N 


The  Fourier  decomposition  of  (3.8)  can  be  done  by  using  the  Fast  Fourier  Transform  (FFT)  of  F, 
which  is  evaluated  numerically  in  the  physical  space. 

The  PSE  equations  (3.7)  can  be  used  to  study  nonlinear  interaction  of  various  modes  (e.g.,  cross¬ 
flow/crossflow,  crossflow/TS,  etc.)  or  one  can  study  the  onset  of  transition  to  turbulence  provided 
appropriate  initial  conditions  are  prescribed.  For  small  disturbances,  F  can  be  neglected  and  one 
obtains  linear  PSE  equations  (after  dropping  the  subscript  11) 


(3.9) 


which  can  be  solved  to  study  the  effect  of  nonparallel  flow  or  that  of  initial  conditions.  If  nonparallel 
effect  is  ignored,  then  (3.9)  essentially  reduces  to  the  Orr-Sommerfeld  equation. 

The  streamwise  wavenumber  in  (3.6b)  needs  to  be  determined  in  order  to  solve  the  equations  by 
a  marching  scheme.  This  procedure  is  given  in  Chang  et  al.  (1991).  Here  we  briefly  describe  it  for 
the  linear  equation  (3.9).  In  this  case,  the  evolution  of  the  shape  function  is  monitored  during  the 
process  of  marching  and  the  wavenumber  is  updated  by  local  iterations  at  a  given  x  according  to  the 
change  in  'P.  At  a  given  location  jci,  let  the  streamwise  wavenumber  be  given  by  Ui  and  then  express 
0as 


^{x,y,z,t)  =  H'(x,y)e 


^^<ndi+fk-ot 


(3.10) 


The  change  of  the  shape  function  V  can  be  approximated  by  the  following  Taylor  series  expansion 
truncated  to  the  first  order 

!P(x.y)=*Pi+^{x-Xi)+... 

OX 

where  'Pi  is  the  shape  function  at  x  =  xi.  To  an  accuracy  of  0(x  -  xi),  the  above  equation  can  be 
further  expressed  as 

•P(x,y)  =  'PiJ..^?"‘'^  (3.11) 

Substituting  (3.11)  into  (3.10),  we  have  the  "effective"  wavenumber  in  the  vicinity  of  Xi  given  by 
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(3.12) 


a  =  ai 


.  1  d'Fi 
‘•Pi  dt 


The  real  part  of  this  effective  wavenumber  represents  the  phase  change  of  the  disturbance  while  the 
imaginary  part  gives  the  growth  rate.  A  disturbance  is  unstable  if  the  imaginary  part  is  less  than 
zero.  Since  the  shape  function  vector  •Pi  depends  upon  y  and  contains  four  dependent  variables  (u,v, 
etc.),  the  value  of  a  computed  by  (3.12)  will  be  a  function  of  the  y  coordinate  and  the  selected 
dependent  variable.  One  can,  for  example,  use  the  shape  function  u  and  the  y  location  where  u 
reaches  its  local  maximum  to  update  the  wavenumber  at  any  given  x  station  as  the  disturbance 
evolves  downstream.  An  alternative  which  is  used  here  is  to  consider  the  following  integral 
condition. 


(3.13) 


which  removes  the  dependence  of  a  on  y.  If  9  is  a  particular  component  of  'Fi  then  the  dependence  of 
a  on  •Pi  is  retained  in  (3.13).  For  three-dimensional  boundary  layers  we  choose  9  to  be  a  vector  with 
components  ( u,v,w).  Equation  (3.13)  is  used  in  tibe  iterative  solution  of  (3.9)  until  the  second  term  in 
(3.13)  vanishes  to  a  prescribed  tolerance.  An  additional  condition  (Chang  et  al.  1991,  Malik  &  Li 
1993)  needs  to  be  satisfied  in  order  to  obtain  solution  of  (3.7)  by  the  space  marching  approach. 

Numerical  solution  of  the  parabolized  stability  equations  requires  discretization  in  both  x  suid  y 
directions.  We  discretize  the  streamwise  derivative  by  a  backward  Euler  step  and  wall-normal 
derivatives  by  fourth-order  accurate  compact  differences  (see,  Mahk,  Chuang  &  Hussaini  (1982)). 
Homogeneous  boundary  conditions  at  the  wall  and  in  the  free  stream  are  imposed.  The  initial 
conditions  are  obtained  by  a  local  approximation  to  (3.9)  and  by  solving  the  associated  eigenvalue 
problem.  Since  the  wave  information  is  absorbed  in  the  wavenumber  a  (3.6b),  one  needs  to  use  a  few 
marching  steps  per  wavelength  to  obtain  an  accurate  solution  of  the  wave  evolution.  Calculations  for 
two-dimensional  boundary  layers  show  that  PSE  results  with  only  3  steps  per  wave  length  agree 
quite  well  with  very  acciuate  Navier-Stokes  computations  using  60  grid  points  per  wavelength  (see, 
Joslin,  Streett  and  Chang  1992). 


4.  Linear  and  Nonlinear  Stability  Analysis  and  Wave  Interaction 

4. 1  Quasi-parallel  linear  stability 

In  order  to  determine  the  relevant  physical  parameter  space,  it  is  appropriate  to  first  give  some 
results  from  quasi-parallel  linear  stability  theory.  We  consider  two  cases:  R  -  250  and  500.  Hail, 
Malik  &  Poll  (1984)  foimd  that  the  attachment-line  boundaiy-layer  (x  =  0)  is  stable  to  infinitesimal 
disturbances  up  to  i2  =  583.1.  Thus,  for  the  two  cases  considered  here,  the  attachment-line 


56 


boundary-layer  is  stable.  In  the  present  study,  we  are  interested  in  the  crossflow  disturbances  which 
will  become  unstable  away  from  the  attachment-line  (x  »  1).  Figure  1  shows  the  mean  velocity 
profiles  in  directions  tangential  and  across  the  inviscid  stream  at  an  R  of  500  and  R  =  500.  The 
velocity  profiles  Ut  and  Uc  are  defined  as 

Ui=u  cos  6 +  w  sin  6 
Uc  =usin6-wcos9 

where  the  streamline  angle  6  is  defined  in  (2.6).  These  velocity  profiles  have  been  scaled  with 
spanwise  inviscid  velocity  Woo-  It  is  clear  from  (2.6)  that  6  decreases  with  i?  or  x  as  the  flow  turns 
away  firom  the  attachment-line  towards  the  free-stream  direction.  This  is  depicted  in  figure  2  where 
the  angle  6,  along  with  the  crossflow  Reynolds  number  (defined  below),  is  plotted  for  both  R  = 
250  and  500. 

Crossflow  instability  is  associated  with  the  inflectional  velocity  profile  Uc  which,  for  swept 
wings,  is  positive  towards  the  center  of  curvature  of  the  streamline.  The  flow  becomes  unstable 
when  crossflow  Reynolds  number  ^  40  where  is  defined  by 


where  Ug  is  the  maximum  value  of  the  crossflow  velocity  Uc  and  d  ^  is  the  thickness  where  the  cross- 

flow  velocity  has  dropped  to  10%  of  Ug.  Distribution  of  the  crossflow  Reynolds  number  R^^  is  given 

for  the  two  cases  in  figure  2.  The  value  of  R^^  exceeds  about  50  at  R  =  200  and,  hence,  the  instability 

will  onset  at  R  <  200  for  both  cases.  The  maximum  value  of  R.,  is  about  150  for  R  =  250  and  about 

'a 

270  for  R  =  500.  In  swept-wing  flows,  transition  usually  occurs  where  R*^  becomes  of  (X200). 

Figure  3  presents  results  for  integrated  growth, 

in{AI  A^)  =  \\pdR  (4.2) 

using  the  quasi-parallel  growth  rate  Op  =  -a,-.  Calculations  are  performed  for  stationary  as  well  as 
traveling  disturbances  with  frequency  R  =  .75x  10~*  (where  F  =  2jzvf  /  W^,  f  being  the  frequency  in 
hertz)  at  both  R .  These  calculations  are  performed  for  a  fixed  spanwise  wavenumber  of  0.4  which  is 
close  to,  but  no  quite  (see  figure  4  below),  the  most  amplified  wave  number  for  the  flow  under  study. 
It  is  clear  that  traveling  disturbances  amplify  more  than  the  stationary  disturbances  according  to 
linear  theory.  However,  stationary  disturbances  are  found  to  dominate  when  experiments  are 
performed  in  low-disturbance  wind  tunnels.  This  is  due  to  the  lower  initial  amplitude  of  traveling 
modes  (see  the  work  of  Choudhari  &  Streett  (1990)  on  the  receptivity  of  stationary  and  traveling 
disturbances). 
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The  variation  of  spatial  growth  rate  with  wavenumber  p  is  given  in  figure  4  for  the  two  frequen¬ 
cies.  There  are  two  curves  associated  with  F  =  .75xl0~*,  one  with  positive  P  and  another  with 
negative  P,  the  latter  with  smaller  growth  rates.  The  stationary  vortex  and  P  >  0  traveling 
disturbance  has  peak  growth  rate  at  ^  .35  as  shown  in  the  figure  for  R  =  300.  In  the  downstream, 

the  peak  shifts  to  higher  wavenumbers  and  lies,  for  example,  atp  ^  .45,  R  =  650.  The  two  families  of 
unstable  traveling  disturbances  are  further  shown  in  figure  5  where  the  growth  rate  of  the  most 
amplified  (among  various  wave  orientations)  disturbance  is  plotted  as  a  function  of  frequency.  The 
family  with  high  growth  rates  has  its  wave  vector  oriented  at  positive  angles  with  respect  to  the 
invisdd  flow  streamline  (angles  measured  from  the  convex  side)  while  the  family  with  lower  growth 
rates  has  its  wave  vector  oriented  at  negative  angles.  The  relative  sense  of  the  two  modes  depends 
upon  the  direction  of  the  crossflow  with  the  more  amplified  mode  always  oriented  opposite  to  the 
crossflow  direction.  In  both  cases,  the  direction  of  the  group  velocity  lies  at  small  angles  to  the 
invisdd  streamline  direction.  Thus,  the  disturbance  energy  propagates  downstream  for  both  modes 
as  also  noted  by  Mack  (1985).  The  traveling  mode  with  lower  growth  rate  may  be  important  in  the 
nonlinear  stage  and  its  interaction  with  the  more  amplified  traveling  mode  may  also  induce 
stationary  crossflow  vortices  when  other  stimuli,  e.g.,  wall  roughness,  are  absent.  Furthermore, 
these  two  traveling  modes  along  with  stationary  vortex  mode  constitute  a  possible  resonant  triad 
which  may  be  relevant  in  the  transition  process. 

4.2  Nonparallel  effects 

We  now  compare  the  quasi-parallel  growth  rate  results  with  those  obtained  by  solving  linear 
PSE  equations  (Eq.  (3.9)).  Figure  6(a)  shows  the  results  for  R  =  250  for  stationary  vortices  while 
figure  6(b)  shows  the  results  for  a  frequency  of  F  =  .75x  10“*.  In  case  of  PSE,  different  growth  rate 
results  are  obtained  for  u,  v  and  w  components  of  velodty.  At  low  R  (between  200  and  400)  there  is 
considerable  difference  between  these  growth  rates  with  u  growth  higher  than  v  and  w  but  the 
latter  two  approach  the  same  value  at  higher  Reynolds  numbers.  Figure  7  shows  the  growth  rate 
results  for  R  =  500.  In  this  case  the  qualitative  trends  are  the  same  but  there  is  less  difference 
between  the  three  growth  rates.  The  quasi-parallel  growth  rate  is,  in  general,  close  to  the  growth 
rate  based  upon  w  component,  except  at  lower  Reynolds  numbers  where  it  lies  somewhere  in 
between  the  three  growth  rates.  Thus,  one  can  not  make  a  strong  statement  about  nonparallel 
effects  except  that  they  are  more  pronounced  at  lower  R  and  that  they  are  destabilizing  if  measured 
by  the  chordwise  velocity  component.  The  growth  rate  can  also  be  defined  based  upon  the  total 
disturbance  energy  which  accounts  for  all  the  velocity  components  and  growth  rates  based  upon  this 
definition  suggests  that  nonparallel  effect  is  usually  destabilizing,  but  there  may  be  some  exceptions. 
Spalart  (1989)  pointed  out  that  the  growth  rates  from  his  simulation  were  very  close  to  the  quasi- 
paraUel  results  and  that  the  agreement  was  better  at  lower  Reynolds  numbers  (R)  than  at  higher 
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Reynolds  numbers.  Figure  7(a)  shows  that  this  is  true  for  the  case  he  ran  ( i?  =  500,  F  =  0),  but  it  is 
not  a  general  statement  as  is  evident  from  the  comparison  of  Figs.  6  and  7. 

Our  results  for  P  =  A  and  R  =  500  are  compared  with  linearized  Navier-Stokes  computation  of 
Streett  (1993)  in  figure  8.  Streett  performed  spatial  simulations  and  solved  the  full  linearized  sys¬ 
tem  where  the  disturbances  were  introduced  by  spanwise  periodic  steady  suction  and  blowing.  After 
the  initial  transients  die  out,  the  agreement  between  the  two  calculations  is  excellent  and  it  remains 
so  for  a  large  chordwise  extent.  Giood  agreement  was  also  found  with  the  results  of  Spalart  (1989) 
(Malik  &  Li  (1992)).  The  agreement  with  full  Navier-Stokes  solution  shows  that  PSE  approximation 
introduces  negligible  error  in  our  study  of  the  crossflow  vortices  as  disturbance  growth  rate  is  a 
sensitive  quantity  and  any  error  would  have  shown  up  in  growth  rate  results. 


4.3  Nonlinear  development  of  stationary  crossflow  vortices 
Navier-Stokes  simulations  by  Malik  (1986)  for  rotating-disk  flow  and  by  Meyer  &  Kleiser  (1988) 
for  a  Falkner-Skan-Cooke  boundary  layer  showed  nonlinear  saturation  of  crossflow  vortices.  Both 
these  calculations  employed  temporal  approach  and,  therefore,  ignored  nonparallel  effects.  Here  we 
present  spatial  nonlinear  calculations  for  R  =  500  using  PSE.  Initial  conditions  for  the  stationar> 
vortex  with  ^  s  .4  were  prescribed  at  R  =  186.  It  was  assumed  that  the  vortex  shape  is  given  by  the 
linear  eigenfunction  at  that  Reynolds  number  and  that  the  maximum  disturbance  amplitude 
(max(«* +u)*)^^*)  is  .001  Figure  9  gives  the  computed  total  perturbation  wall  vorticity 
distribution 


which  shows  streamwise  striations  starting  at  i?  =  350.  The  green  color  indicates  negative  values 
while  the  red  indicates  positive.  The  perturbation  wall  vorticity  values  are  very  small  initially  and 
the  signal  becomes  noticeable  (strong)  only  atiZ  of  about  420.  As  we  will  show  later,  the  disturbance 
amplitude  at  this  location  has  already  reached  about  4  percent.  Hence,  when  these  vortices  are 
observed  in  a  flow  visualization  experiment  it  is  sdmost  certain  that  they  have  entered  the  nonlinear 
stage  with  growth  rates  somewhat  smaller  than  that  given  by  the  linear  theory.  These  striations  are 
evident  in  almost  all  crossflow  experiments  (Gray  (1952),  Gregory,  Stuart  &  Walker  (1955),  Poll 
(1985),  Saric,  Daganhart  &  Mousseux  (1989))  and  result  due  to  variation  in  the  wall  shear  caused  by 
stationary  vortices.  These  vortices  make  a  small  angle  (4-5°)  with  respect  to  the  inviscid  free  stream. 

Figure  10  shows  the  contours  of  u  velocity  in  y-z  plane  at  various  Re3molds  numbers  (R  =  400, 
500,  600  and  650).  Two  spanwise  wavelengths  are  shown  and  the  y  coordinate  has  been  stretched  for 
clarity.  Crossflow  vortices  appear  to  result  in  a  half-mushroom-like  structure  which  is  shown 
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exaggerated  in  the  flgure.  The  actual  structure  is  much  more  flat  as  bhown  in  figure  11  drawn  to 
scale.  Initially  the  boundary-layer  thickness  is  constant  in  z;  however,  as  the  crossflow  instability 
rolls  up  into  vortices,  there  appear  regions  of  low  and  high  velocity  and,  therefore,  the  boundary- 
layer  thickness  varies  considerably  in  the  span  as,  for  example,  seen  for  R  =  600.  In  this  case  the 
variation  is  as  much  as  by  a  factor  of  about  4. 

There  is  a  region  near  z  =  10  and  25  (for  R  =  6(X))  where  the  fluid  is  pushed  towards  the  wall 
while  it  is  pushed  away  from  the  wall  near  z  =  5  and  20.  It  is  these  low  velocity  regions  at  z  =  5  and 
20  where  oil  accumulates  in  a  flow  visualization  experiment  resulting  in  wall  streaks  such  as  those 
shown  in  figure  9.  The  half-mushroom  structure  observed  in  figure  10  is  the  result  of  the  asymmetry 
induced  by  the  crosswind.  In  two-dimensional  flow  over  a  concave  wall  which  is  subject  to 
centrifugal  instability,  a  full  mushroom  structure  appears  as  experimentally  observed  by  Swearingen 
&  Blackwelder  (1987)  and  Peerhossaini  &  Wesfreid  (1988). 

Figure  12  shows  a  velocity  vector  plot  in  the  y-z  plane  at  fixed  R  of  650.  The  velocities  have 
been  projected  onto  a  cross-section  normal  to  the  vortex  axis.  An  insight  into  the  crossflow  vortex 
structure  may  be  achieved  by  releasing  die  particles  at  some  location  within  the  flow  field  and  fol¬ 
lowing  their  paths  as  they  are  carried  through  the  fluid  in  the  y-z  plane.  Two  particles  are  injected 
at  about  z  =  22  but  one  is  released  very  near  the  wall  while  the  other  is  released  at  y  =  1.7.  The 
latter  particle  rolls  into  a  big  vortex  centered  at  about  y  =  2.5  and  z  =  12.  This  is  the  primary  cross- 
flow  vortex.  There  is  a  second  tiny  vortex  near  the  wall  centered  at  about  y  =  1  and  z  =  8  to  which 
the  particle  released  near  the  wall  is  attracted.  This  second  vortex  which  was  much  weaker  atR- 
600  has  also  been  observed  by  R.-S.  Lin  (private  communications,  1992)  in  his  Navier-Stokes 
simulations  of  the  crossflow  vortex  on  a  swept-wing.  It  should  be  stressed,  however,  that  the  actual 
flow  is  fuUy  three-dimensional  and  varies  along  x.  Hence,  these  particle  traces  do  not  depict  the 
three-dimensional  physical  picture  and  have  been  used  merely  to  facilitate  the  visualization  of  the 
crossflow  vortices. 

Contours  in  figure  KKc)  show  a  second  low  velocity  region  near  z  =  15  and  30;  a  hot  wire  located 
at  y  =  1,  for  example,  will  show  two  velocity  defects  per  wavelength  when  traversed  in  the  spanwise 
direction.  This  is  depicted  in  figure  13  which  shows  that  the  second  defect,  caused  by  the  2p  mode 
and  sometimes  referred  to  as  vortex  doubling,  is  much  smaller  than  that  caused  by  the  main  vortex. 
In  our  simulations,  the  2j3  mode  is  excited  through  nonlinear  interaction  and  its  amplitude  remains 
smaller  than  the  primary  mode  with  wavenumber  P  as  shown  in  figure  14  where  the  amplitude 
functions  for  the  stationary  vortex  along  with  its  harmonics  and  mean  flow  correction  are  plotted  at 
R  =  600.  In  a  laboratory  experiment,  the  2P  mode  may  be  excited  via  surface  imperfections  and  the 
relative  amplitude  of  ip  and  2p  modes  may  be  different  from  the  present  case.  The  disturbance 
amplitude  atR  =  600  has  reached  to  about  30  percent  (when  scaled  with  W„)  with  the  maximum 
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meanflow  correction  of  about  15  percent.  In  a  laboratory  experiment,  this  picture  will  be  altered  due 
to  possible  secondary  instabilities  and  interaction  with  traveling  modes. 

Figure  15  shows  the  velocity  profile  along  the  crossflow  vortex  at  4  different  locations  across  it 
for  R  =  500.  The  base  flow  given  by  the  Hiemenz  problem  is  sdso  included.  It  can  be  seen  that  these 
profiles  become  strongly  inflectional  due  to  the  motion  within  the  crossflow  vortex.  Such  profiles 
were  also  noted  in  the  swept-wing  experiments  of  Dagenhart  et  al.  (1989)  and  Muller  &  Bippes 
(1988).  These  inflectional  profiles  as  well  as  the  inflectional  profiles  in  z  (figure  13)  are  subject  to 
inviscid  secondaiy  instabilities  which  are  most  likely  related  to  the  high  frequency  disturbances 
observed  by  Kohama,  Saric  &  Hoos  (1991).  We  will  investigate  this  aspect  of  the  problem  in  a  later 
section.  Here,  we  first  consider  the  interaction  of  travehng  and  stationary  crossflow  disturbances. 

4.4  Statiorjiry  and  traveling  wave  interaction 

It  was  pointed  out  that  the  experiment  of  Miiller  &  Bippes  (1988)  suggests  an  early  nonlinear 
interaction  between  stationary  and  traveling  waves.  We  now  consider  such  interactions  in  the 
swept-Hiemenz  flow.  Our  calculations  are  performed  using  =  A  for  both  the  stationsu^  and 
traveling  (F  =  .75xl0~*)  disturbances.  The  initial  conditions  were  imposed  atF  =  186  and  the  ampli¬ 
tude  of  the  stationaiy  wave  was  the  same  as  in  §  4.3  above,  i.e.,  .1  percent  For  traveling  waves,  two 
different  initial  amplitudes  were  considered:  .01  percent  and  .1  percent.  Results  for  both  the  cases 
are  discussed  below. 

Figure  16  gives  the  results  of  disturbance  energy  of  various  modes  denoted  as  (0,1),  (1,1),  (2,2), 
etc.  Here,  the  first  index  refers  to  frequency  to  and  the  second  index  to  spanwise  wavenumber  p. 
Thus,  mode  (2,2)  is  the  harmonic  with  twice  the  frequency  and  twice  the  wavenumber  of  the 
traveling  mode.  For  comparison,  the  case  of  stationary  vortex  only  is  also  given.  For  the  stationary 
vortex  case,  shown  in  figure  16(a),  the  energy  cascades  into  2p,  3P,  4)3 ...  modes  as  earlier  noted  in  the 
simulations  by  Malik  (1986)  and  Meyer  &  Kleiser  (1988).  The  energy  in  the  mean  flow  correction 
mode  is  of  the  same  order  as  the  2p  mode.  It  is  probable  that  the  essential  features  of  the  nonlinear 
development  of  the  stationary  crossflow  vortex  can  be  captured  by  a  model  which  considers  0,  p  and 
2P  modes. 

The  interacting  case  with  stationary  vortex  amplitude  10  times  higher  than  traveling  is  shown 
in  figure  16(b).  This  case  is  meant  to  simulate  moderately  low  turbulence  conditions  where  wall 
roughness  will  introduce  dominant  instability  (i.e.,  stationary  vortex)  and  the  weak  turbulence  will 
introduce  traveling  disturbances  with  low  amplitude.  On  the  other  hand,  figure  15(a)  can  be  thought 
of  as  the  case  with  ultra-low  turbulence  with  essentially  no  traveling  modes  induced.  In  contrast, 
figure  16(c)  is  the  high-turbulence  case  where  the  initial  amplitude  of  the  traveling  mode  is  equal  to 
the  stationary  mode.  Admittedly,  these  are  all  idealized  cases,  for  in  natural  environment  energy 
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input  is  into  a  broad  band  of  frequencies  and  wavenumbers  which  we  cannot  attempt  to  tackle  in  the 
present  framework. 

Figure  16(b)  shows  that  the  energy  in  the  stationary  as  well  as  the  traveling  mode  saturates  at 
about  R  =  490,  the  energy  in  the  latter  mode  remains  smaller  except  near  the  very  end  at  about  R  = 
600  where  the  two  become  the  same.  This  is  also  where  the  energy  in  mode  (l.-l)  supersedes  the  two 
primary  modes  and  becomes  comparable  to  the  mean  flow  correction  mode.  The  mode  (1,-1)  is 
generated  due  to  interaction  between  (1,1)  and  (0,2)  mode  and,  apparently,  the  combination  of  the 
amplitudes  is  just  about  right  to  yield  a  resonance  between  the  three  modes  as  speculated  in  §  4.1 
above. 

The  situation  changes  when  the  initial  amplitude  of  the  two  waves  become  the  same.  The 
traveling  mode  has  the  higher  energy  all  the  way  and  it  tends  to  suppress  the  growth  of  the 
stationary  modes.  Both  the  primary  modes  saturate  earlier  at  about  R  =  430  as  compared  to  490  in 
figure  16(b).  The  higher  modes  gaining  the  dominant  energy  appear  to  be  (2,2),  (3,3),  (4,4)...  modes 
in  this  case.  The  suppression  (also  compare  Figs.  17(a)  and  (c)  below)  of  the  stationary  vortices  by 
traveling  modes  is  supported  by  the  observation  made  in  the  DLR  experiment. 

The  evolution  of  the  maximum  (in  y)  disturbance  amplitude  for  the  three  cases  is  given  in  figure 
17  on  a  natural  log  scale.  Amplitudes  of  all  the  velocity  components  are  given.  Initially  the  spanwise 
velocity  w  is  higher  than  the  chordwise  velocity  u .  Later  the  magnitude  of  the  two  switches  as  the 
inviscid  streamline  angle  decreases  (note  that  6  =  45°  when  R  s  500).  The  magnitude  of  the  normal 
velocity  v  is  much  lower  than  u  and  w  for  both  waves  at  all  R.  From  figure  17(a)  for  stationary 
vortices  alone,  it  is  clear  that  the  nonlinear  N  factors  {in  A/ Aq)  at  R  =  650  are  7  and  5  for  u  and  w. 
These  are  to  be  contrasted  with  the  value  of  about  9  given  by  quasi-parallel  linear  calculations  in 
figure  3. 

The  growth  rates  of  the  stationary  and  traveling  waves  for  the  above  three  cases,  along  with  an 
additional  case  of  traveling  mode  only  (initial  amplitude  of  0.1  percent),  are  given  in  figure  18.  For 
comparison,  the  growth  rate  from  linear  PSE  calculations  are  also  given.  This  plot  more  clearly 
shows  the  behavior  of  the  two  modes  discussed  with  reference  to  figure  16.  The  growth  rate  of  the 
stationary  vortex  (curve  2)  begins  to  depart  firom  linear  theory  at  about  R  =  420.  At  this  location  the 
distiu-bance  amplitude  is  only  about  4  percent  At  i?  =  450,  the  growth  rate  is  lower  than  the  linear 
theory  result  by  about  9  percent,  but  it  begins  to  decrease  rapidly  beyond  that.  The  resulls  are 
similar  for  the  traveling  mode  alone  (curve  2)  with  initial  amplitude  of  .1  percent  Since  the  trav¬ 
eling  mode  amplifies  more  rapidly,  it  reaches  saturation  earlier  and  its  growth  rate  begins  to  depart 
from  the  linear  theory  results  at  i?  =  330.  For  the  wave  interaction  case  with  At  =  .01  percent,  the 
growth  rate  of  the  two  waves  begins  to  depart  fi’om  the  linear  theory  result  at  about  R  =  390.  At 
R  =  450  the  two  growth  rates  differ  from  the  linear  results  by  about  18  percent .  Subsequently,  the 
two  growth  rates  drop  sharply  and  at  12  =  500,  the  stationary  and  traveling  disturbance  growth  rates 


62 


are  lower  by  about  70  and  76  percents  with  respect  to  their  linear  growth  rates.  Hence,  the  results 
indicate  that  even  for  the  case  with  smaller  initial  traveling  disturbance  amplitude  there  is  some 
interaction  well  before  R  of  about  500.  This  interaction  becomes  stronger  when  the  initial 
amplitudes  of  the  two  waves  are  the  same  (A,  =  A|  =  .1  percent).  In  this  case  the  growth  rates  begin 
to  depart  from  the  linear  theory  results  at  i2  =  330  and  by  410  the  growth  rates  have  dropped  by  68 
percent  for  the  stationary  vortex  and  by  39  percent  for  the  traveling  mode.  A  close  examination  of 
the  results  show  that  when  A(  =  .1%,  there  is  no  direct  effect  of  the  stationary  disturbance  on  the 
traveling  wave  (curves  3  and  5  collapse)  but  the  growth  of  the  stationary  vortex  is  greatly  suppressed 
due  to  the  presence  of  higher  amplitude  traveling  disturbance.  It  is  clear  that  the  two  modes  do 
interact  depending  upon  the  initial  amplitude,  as  also  inferred  by  Bippes  (1991)  from  his 
experiments.  The  nonlinear  growth  rate  behavior  at  large  R  indicates  that  the  two  primary  modes 
reach  a  quasi-equilibrium  state  where  the  growth  rate  begins  to  oscillate  around  a  small  value. 

In  order  to  shed  some  more  light  on  the  stationary/traveling  mode  interaction,  we  consider  the 
case  with  A,  =  .1  percent.  At  =  .01  percent,  and  plot  in  the  y-z  plane  at  four  different  Reynolds 
numbers  (R  =  431,  500  and  600).  Here  u^is  defined  as 
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The  results  are  shown  in  figure  19  for  the  three  locations.  Due  to  nonlinearity  and  interaction 
with  the  stationary  mode,  traveling  disturbances  are  modulated  in  the  spanwise  direction.  An 
important  observation  is  that  the  peak  rms  perturbation  is  near  the  wall  at  y  ~  1  with  a  second 
maximum  (but  with  much  lower  amplitude)  further  away  from  the  wall  (see  figure  19(b)).  A 
comparison  of  Figs.  19(b)  and  10(b),  at  R  =  500,  shows  that  the  peak  occurs  in  the  spanwise 
region  where  low-velocity  fluid  is  pushed  away  from  the  wall.  At  higher  Reynolds  number  (I?  =  600), 
there  are  two  peaks  in  near  the  wall,  apparently  associated  with  the  emergence  of  2)3  harmonic 
of  the  stationsuy  mode.  Michel,  Amed  &  Juillen  (1985)  also  noted  two  maxima  in  the  root  mean- 
square  value  of  the  streamwise  velocity  within  a  spanwise  wavelength  in  the  ONERA/CERT  swept- 
wing  experiment.  The  magnitude  of  the  maximum  value  was  found  to  be  up  to  about  20  percent  of 
the  resultant  inviscid  velocity.  They  also  found  that  most  of  the  turbulence  energy  is  contained  in 
the  firequency  range  which  is  unstable  according  to  the  linear  stability  analysis. 

Figure  20  is  a  plot  of  the  stationary  as  well  as  rms  velocity  signal  («  component)  aty  =  1.048. 
Modulation  of  the  traveling  disturbances  due  to  the  presence  of  stationary  vortex  is  evident.  The 
peak  is  in  the  region  where  a  velocity  defect  appears  in  the  stationary  signal  and  the  minimum 
in  occurs  where  there  is  a  velocity  access.  However,  there  is  a  phase  shift;  of  about  fl/4  between 
the  maximum  in  and  the  minimum  stationary  velocity,  as  evident  from  results  for  R=  500.  This 
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phase  shift  decreases  at  higher  Reynolds  numbers.  Muller  &  Bippes  (1988)  reported  experimental 
results  qualitatively  similar  to  those  in  figure  20. 

We  plot  the  variation  of  with  Reynolds  number  aty  =  1.048  in  figure  21  which  shows  the 
maximum  and  minimum  as  well  as  along  the  path  where  stationary  u  velocity  has  a 
maximum  and  a  minimum.  The  variation  of  the  stationary  u  velocity  with  Reynolds  number  is  also 
shown  in  the  figure.  The  figure  clearly  shows  that  in  high  (stationary)  velocity  region  the 
component  gets  satiu'ated  but  it  increases  to  higher  amplitudes  in  the  low  (stationaiy)  velocity 
region.  At  R  =  550,  the  maximum  reaches  about  20  percent  in  the  low-velocity  region.  Such 
high  levels,  although  dependent  upon  the  initial  disturbance  amplitude,  are  not  unexpected  in  view 
of  the  experimental  evidence  provided  by  Michel,  Amal  &  Juillen  (1985).  Poll  (1985)  also  observed  a 
traveling  disturbance  with  frequency  close  to  the  most  amplified  disturbance  given  by  linear  theory. 
He  further  noted  that  close  to  the  surface  the  amplitude  of  these  disturbances  can  exceed  20  percent 
of  the  local  mean-flow  velocity.  However,  the  rms  amplitudes  measured  by  Dagenhart  et  al.  (1989) 
are  much  lower  which  suggests  that  in  their  experiment  the  initial  amplitude  of  traveling  modes 
relative  to  the  stationary  mode  was  much  lower  than  used  here.  Choudhari  (1993)  estimates  that 
the  initial  amplitude  of  the  traveling  mode  could  be  up  to  two  orders  of  magnitude  lower  than  the 
stationaiyr  modes  for  the  receptivity  mechanism  considered  in  his  study.  This  may  possibly  be  the 
case  in  the  experiment  of  Dagenhart  et  al. 

4.5  Effect  of  nonlinear  disturbances  on  skin-friction 

Figure  22  gives  the  chord  wise  (C^)  and  span  wise  skin-friction  coefficients  for  all  three 
cases.  From  figure  9  we  know  that  skin-friction  varies  in  the  spanwise  direction.  However,  figure  22 
gives  the  spanwise  averaged  value,  i.e.,  only  the  contribution  from  mean  flow  distortion  is  consid¬ 
ered.  Since  w  is  independent  of  R,  lamin^u'  spanwise  skin-friction  remains  constant.  Similarly, 
since  u  increases  linearly  with  R,  so  does  the  chordwise  skin-firiction  when  scaled  with  W^.  At  some 
location  both  and  begin  to  depart  from  their  respective  laminar  values.  For  case  (a)  and  (b) 
of  figure  16,  this  location  is  at  about  R  =  450  and  from  thereon  it  rises  significantly.  The  skin  friction 
rise  from  stationary  vortex  alone  is  about  19  percent  for  and  58  percent  for  C/„,  at  R  =  600.  The 
skin-firiction  for  case  (b)  is  slightly  higher  in  the  beginning  but  later  on  it  drops  below  case  (a). 
Computations  for  case  (a)  were  made  using  N  =  2,9  and  16  in  (3.5).  While  there  was  considerable 
difference  in  the  skin-friction  distribution  for  iV  =  2  and  9,  essentially  no  difference  was  found 
between  the  two  higher  resolution  cases.  For  cases  (b)  and  (,c)M  =  N  =  9  was  used  in  (3.5). 

In  case  (c),  with  higher  initial  amplitude  of  the  traveling  mode,  skin  friction  begins  to  rise  much 
earlier  at  about  R  =  375  as  would  be  expected  from  the  comparison  of  figiu*e  16(b)  and  (c)  which 
shows  that  the  mean  flow  distortion  is  higher  in  the  latter  case.  Hence,  a  stronger  interaction  of  the 
traveling  and  stationary  modes  leads  to  higher  skin-friction  coefficient.  Our  results  indicate  that  the 
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angle  between  wall  shear  and  inviscid  free  stream  decreases  as  the  disturbed  flow  enters  highly 
nonlinear  stage  in  the  three-dimensional  boundary  layer. 

5.  High  Frequency  Secondary  Inatability 

A  hot-wire  placed  in  a  three-dimensional  boundary  layer,  subject  to  crossflow  instability,  sees 
two  types  of  unsteady  disturbances.  First,  it  captures  an  unsteady  signal  with  a  peak  at  a  frequency 
which  coincides  with  the  most  amplifled  frequency  given  by  the  linear  stability  theory.  In  the 
present  case,  e.g.,  /i  =  1.2  x  10“*  at  R  =  300  (see,  figure  b).  The  dimensional  value  of  this  frequency 
depends  upon  the  flow  parameters  (unit  Reynolds  number,  sweep,  etc.).  Poll  (1985),  for  his  swept 
cylinder  experiment,  found  /i  to  be  1500  Hz  for  the  chord  Reynolds  number  of  1.2x10®  and  sweep 
angle  of  63°.  In  Kohama,  Saric  &  Hoos  (1991),  was  close  to  180  Hz  and  in  Michel,  Amal  &  Juillen 
(1985)  fi  <  200.  In  these  experiments,  the  hot-wire  also  captures  a  second  frequency  /“j  which  is  an 
order  of  magnitude  larger  than  For  example,  was  17500  Hz,  3500  Hz  and  1000  Hz  in  the 
experiments  of  PoU  (1985),  Kohama,  Saric  &  Hoos  (1991)  and  Michel,  Arnal  &  Juillen  (1985), 
respectively.  In  this  section,  we  investigate  this  high  frequency  instability  in  the  present  three- 
dimensional  boundary  layer.  The  problem  is  modeled  here  as  the  secondaiy  instability  of  the  new 
mean  flow  which  is  set  up  by  the  presence  of  a  large-amplitude  stationary  crossflow  vortex. 

We  perform  secondary  instability  analysis  locally,  i.e.,  at  a  fixed  Reynolds  number  and  perform 
temporal  stability  analysis.  In  order  to  perform  this  analysis,  we  rotate  the  x-z  coordinates  to  a  new 
system,  X2,  so  that  the  X2  coordinate  aligns  with  the  crossflow  vortex.  At  a  streamwise  location 
designated  by  the  Reynolds  number  R,  we  ignore  the  curvature  of  the  vortex  and  use  the  quasi¬ 
parallel  approximation  (we  will  provide  a  posteriori  justification  for  these  assumptions  later)  which 
allows  us  to  consider  a  harmonic  disturbance  of  the  type 

where  a2  and  <02  are  the  wavenumber  and  frequency  of  the  secondary  disturbance  and  y2  =  y-  Here, 
since  we  use  temporal  stability  concept,  is  real  and  (Wg  is  complex.  If  Og,  >0  (a)2x  =Iinag(ffl2)). 
then  the  secondary  instability  is  present.  Temporal  stability  approach  has  earlier  been  used  by 
Herbert  (1983)  for  secondary  instability  of  TS  waves  and  by  Hall  and  Horseman  (1991)  for  secondary 
instability  of  (Jdrtler  vortices.  This  approach  can,  at  least,  provide  a  qualitative  picture  of  the 
secondary  instability  phenomenon. 

We  superimpose  (5.1)  on  the  meanflow  computed  in  §  4.3  above,  i.e.,  the  meanflow  constitutes 
the  three-dimensional  boundary  layer  as  modulated  by  the  presence  of  a  nonlinear  stationary 
crossflow  vortex  with  initial  amplitude  of  .1  percent.  This  meanflow,  when  represented  in  (X2,y2»^2) 
coordinate  system,  is  a  strong  function  of  yg  and  Z2  but  a  weak  function  of  x^-  Substituting  the 
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meanflow  and  the  disturbance  wave  (5.1)  in  incompressible  Navier-Stokes  equations,  we  obtain  the 
following  linearized  equations 


ia2f/2+a|/i?  +  ^  U2+^U2+^M/2+i«2P2+'^2^  +  W'2^--i  +  =  i032U2  (5.2) 

COC2  J  <^2  ^2  ^2  ^  |_  ^2  ^2  J 

3V2  r.  rr  ,  2/n  <^2!  <^2  xr  <^2  irr  <^2  4^2  1  r<?^y2  1 

^U2+  ia2U2+a2lR  +  -^  v^+-^Wi  +  V2-^  +  W^-^  +  -^-—  -^  +  —^  =10)202  (5.3) 

OT2  [  ^2  J  <®2  ^2  *2  <?y2  ^  L  ^2  *2 

dWo  dWo  rr  .  2  in  ^2  xr  ^Wn  dWn  dp,,  1  d^Wo  d^Wo  .... 

^02+3-^02+  ia2U2+a^lR  +  ^  (02+V2-^+W2^  +  ^-—  ^-#  +  ^-2^  =10)2^2  v5.4) 

®*2  L  “2  J  <^2  *"2  ®2  L  ^2  J 


.  ^^2  dw^ 

ia2U2  +Tr^  +  ^^  =  0 
£^2  dZ2 


where  U2,  V2,  W2  are  the  mean  velocity  components  in  X2,  y^,  Z2  directions,  respectively,  and 
U2.  02,  W2  are  the  corresponding  disturbance  velocity  components  and  P2  is  the  pressure.  In  (5.2)- 
(5.4)  terms  dU2ldx2,  3^2! dx^,  dW2ldx2  are  small  and  can  be  neglected  as  numerical  experiments 
indicate  that  they  do  not  appreciably  change  the  eigenvalue.  However,  ^2 1^2  is  of  the  same  order 
as  dW2l  dz2  and  thus  V2  must  not  be  set  to  zero.  Dropping  V2  increases  the  growth  rate  by  about  50 
percent.  Equations  (5. 2-5.5)  are  partial  differential  equations  which  are  subject  to  homogeneous 
conditions  at  the  wall  and  free-stream,  i.e., 

Ug  =  Ug  =  0,2  =  0,  ^2  =  0  (5.6a) 

«2  -)  0,  V2-*0,  W2->0  as  y2  (5.6b) 

The  computational  domain  in  22  direction  covers  one  wavelength  of  the  stationary  vortex  and 
periodic  boundary  conditions  are  imposed  in  the  Zj  direction,  i.e.. 


“2(22)  =  “2(22 

(5.7a) 

*^2(22) =*^2(^2 +4) 

(5.7b) 

0)2(22)  =  102(^2 +-^7)  - 

(5.7c) 

where  =  2)r/(af  +  a,  and  being  the  x  and  z  wavenumbers  of  the  stationary  vortex. 

Equations  (5.2-5.5)  along  with  the  boundary  conditions  (5. 6-5. 7)  constitute  an  eigenvalue  prob¬ 
lem  which  we  solve  by  using  a  Chebyshev  collocation  method  in  the  >2  direction  and  a  Fourier  collo¬ 
cation  method  in  the  Z2  direction.  The  physical  domain  yo  €[0,y2  1  is  mapped  on  to  a  computa¬ 

tional  domain  tj  e  [-1, 1]  such  that  the  grid  points  are  clustered  near  the  wall  and  ^2  =  ^21  >  where  y2 
is  the  location  where  the  secondary  structure  is  concentrated.  Since  we  do  not  stagger  the  mesh  in 
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direction,  two  additional  boundary  conditions  are  required  which  we  prescribe  by  evaluating  the 
normal  momentum  equation  at  y2  =  0  >2^  • 

The  eigenvalue  problem  can  be  represented  in  the  form 

(5.8) 

where  B  is  a  diagonal  matrix  and  A  is  a  square  matrix  of  size  -  6).^^  where  Ny  and  N,  are  the 
number  of  collocation  points  in  y2  and  22  directions,  respectively.  The  above  eigenvalue  problem 
(5.8)  is  solved  by  the  QR  method  which  yields  all  the  eigenvalues  of  the  discretized  system.  We  test 
the  accuracy  of  these  eigenvalues  by  using  inverse  Rayleigh  iteration  method.  Among  the  computed 
eigenvalues,  only  a  few  have  CO2.  >  0 .  Here,  we  discuss  only  one  of  these  eigenvalues. 

Computations  were  first  performed  at  R  =  450  where  the  stationary  crossflow  disturbances  had 
gained  an  amplitude  of  about  8  percent  based  upon  the  local  inviscid  velocity  (see,  figime  17(a)).  No 
secondary  instability  was  found,  i.e.,  all  0)2.  <0.  The  analysis  was  then  repeated  at  R  =  500  where 
the  maximum  u  amplitude  is  about  17  percent  and  the  stationary  vortex  is  on  its  way  to  saturation 
(compare  with  figure  18).  At  this  location  the  secondary  instability  is  found  but  the  growth  rate  is 
about  the  same  as  that  of  the  nonlinear  crossflow  vortex.  Finally,  calculations  were  performed  at  R  = 
550  where  figure  18  shows  that  the  stationary  vortex  has  a  growth  rate  which  is  close  to  zero.  The 
local  maximum  amplitude  of  the  stationary  vortex  is  about  22  percent.  Secondary  instability  results 
for  this  Reynolds  number  are  discussed  next. 

The  frequency  (0)2,)  of  the  secondary  instability  and  temporal  growth  rate  (0)21)  are  plotted  in 
figure  23.  The  peak  growth  rate  of  C02i  ~  .02  occurs  at  a2  of  about  .6.  At  this  location  (02r  is  about 
.75  which  amounts  to  an  F2  of  1.5x10“^.  We  noted  earlier  that  the  most  amplified  traveling 
crossflow  disturbance  has  a  frequency  of  about  Fi  =  1.2x  10~*;  hence,  F2  is  an  order  of  magnitude 
higher  than  F^  which  is  in  agreement  with  the  experiments  mentioned  above.  The  wavenumber  of 
the  stationaiy  disturbance  along  Z2  coordinate  is  about  .5.  Given  that  02  is  .6,  the  angle  of  the 
secondary  structure  is  about  50°  with  respect  to  the  crossflow  vortex.  The  secondary  instability 
convects  along  the  stationary  vortex  with  the  phase  velocity  of  about  1.25  Woo.  The  relatively  high 
amplitude  of  the  stationary  vortex  required  for  secondaiy  instability  is  in  agreement  with  the  results 
of  Balachandar,  Streett  &  Malik  (1992). 

The  above  calculations  were  made  with  Ny  =  41  and  =  8.  Since  8  collocation  points  in  the 
spanwise  direction  may  be  too  few,  we  repeated  some  of  the  calculations  with  =  16  and  Ny  =51. 
These  results  are  also  given  in  figure  23.  We  note  that  although  there  is  some  movement  in  the 
eigenvalues,  the  results  given  in  the  figure  with  the  lower  resolution  are  qualitatively  correct,  at 
least  at  high  a2-  At  lower  values  of  the  wavenumber  a2,  there  appears  to  be  an  intricate  mode 
structure  which  for  its  investigation  would  require  the  development  of  more  efficient  means  of 
computing  eigenvalues  of  very  large  matrices  so  that  these  calculations  can  be  readily  performed. 
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The  quasi-parallel  approximation  used  in  the  above  analysis  can  be  justified  since  our  results 
show  that,  in  the  rotated  coordinate  system,  variations  of  f/2»^2>^2  along  is  much  smaller  than 
the  variation  in  y2,  Z2  and  the  X2-wavelength  of  the  secondary  instability  is  only  about  3  times  the 
boundary-layer  thickness.  In  the  nondimensional  units,  this  wavelength  for  the  most  amplified 
secondary  wave  is  about  10.  The  secondary  instability  analysis  shown  in  Fig.  23  is  performed  at  A  = 
550.  It  is  clear  from  Fig.  9  that  the  curvature  of  the  vortex  in  the  range  of  iZ  =  550  ±  5  can  be 
ignored. 

The  structure  of  the  secondary  instability  is  presented  in  figure  24  where  the  \i^\  eigenfunction 
is  plotted  along  with  the  U2  component  of  the  meanflow  in  I  =  550.  It  is  clear  that 

this  high-frequency  instability  resides  on  top  of  the  stationary  vortex  with  the  maximum  in 
located  near  ^  2.8 .  In  contrast,  the  lower  frequency  traveling  crossflow  disturbance  is 
concentrated  near  the  wall  at  y2  =  1  (see  figure  19).  The  high  frequency  instability  is  inviscid  in 
nature  and  it  can  be  captured  by  dropping  viscous  terms  in  (5.2-5.4);  however,  since  the  basic  flow  is 
three-dimensional  and  varies  with  y2  and  Z2  it  is  not  possible  to  reduce  the  problem  to  a  single 
partial  differential  equation.  However,  numerical  experiments  suggest  that  some  qualitative 
features  of  the  instability  can  be  captured  by  considering  just  the  component  of  the  mean  flow  in 
which  case  a  single  partial  differential  equation  can  be  used  resulting  in  substantial  savings  in 
computer  time. 

The  top  view  of  the  flow  field  that  results  by  superimposing  the  secondary  eigenstructure  (with 
an  amplitude  of  5%)  on  the  meanflow  ((/2,V2,W2)  is  depicted  in  figure  25  where  the  X2-velocity 
component  is  plotted  at  y2  =  2.82.  Two  periods  in  both  and  Z2  directions  are  shown.  The  dark 
patches  in  the  picture  correspond  to  the  corotating  structures  which  move  in  the  x^  direction.  A  hot 
wire  placed  near  the  boundary-layer  edge  will  detect  this  high-frequency  disturbance  but  if  the  hot¬ 
wire  is  located  in  the  region  between  ^2  =  5  and  10,  for  example,  this  instability  will  not  be  captured. 
Therefore,  extreme  care  is  needed  in  order  to  detect  this  secondary  structure  in  an  experiment. 

6.  Conclusions 

We  have  investigated  crossflow  instability  in  a  model  three-dimensional  boundary  layer  which 
has  an  exact  solution  to  the  incompressible  Navier-Stokes  equations.  This  consists  of  the  swept- 
Hiemenz  flow  which  forms  near  an  attachment-line.  This  flow  is  subject  to  Tollmien-Schlichting 
instability  for  small  x,  where  x  is  the  chordwise  distance,  provided  the  spanwise  Re3aiolds  number  R 
>  583.1.  However,  this  boundary-layer  becomes  unstable  to  crossflow  instability  for  x  »  1  even  for 
R  <  583.1.  Here,  we  have  considered  R  =  250  and  500  for  the  linear  stability  and  R  =  500  for  the 
nonlinear  case.  Both  the  linear  and  nonlinear  stability  as  well  as  the  wave-interaction  in  this  three- 
dimensional  boundary  layer  is  studied  using  parabolized  stability  equations  (PSE).  We  also  study 


68 


secondary  instability  in  this  boundary  layer.  We  find  that  the  various  features  of  the  swept-wing 
boundary-layer  transition  are  captured  in  the  study  of  this  model  boundary  layer. 

Our  linear  results  show  that  nonparallel  effects  are  destabilizing  for  crossflow  disturbances. 
However,  the  magnitude  of  the  effect  depends  upon  R  (more  destabilization  at  lower  R ).  Growth 
rate  of  stationary  crossflow  vortices  computed  from  linear  PSE  are  in  agreement  with  the  results 
obtained  from  Navier-Stokes  simulation. 

Nonlinear  development  of  stationary  crossflow  vortex  is  also  investigated  for  an  initial 
amplitude  of  0.1  percent.  The  computed  wall  vorticity  distribution  shows  the  familiar  streamwise 
streaks,  in  agreement  with  the  surface  oil-flow  visualizations  in  swept-wing  experiments.  Other 
features  of  the  stationary  vortex  development  observed  in  the  experiments  (half-mushroom 
structure,  highly  inflected  velocity  profiles,  vortex  doubling,  etc.)  are  also  captured  in  our  nonlinear 
PSE  computations. 

Nonlinear  interaction  of  stationary  and  traveling  crossflow  modes  is  also  studied.  When  the 
initial  amplitude  of  stationary  vortex  is  large  as  compared  to  the  traveling  mode,  the  stationary 
vortex  dominates  most  of  the  downstream  development.  Eventually,  however,  the  traveling  mode 
becomes  of  the  same  order  as  stationary  mode.  Interaction  of  the  traveling  mode  with  the  harmonic 
of  the  stationary  mode  gives  rise  to  another  traveling  mode  with  same  frequency  but  negative 
spanwise  wavenumber.  Apparently,  a  triad  resonance  is  set  up  at  this  stage.  The  situation  changes 
when  the  initial  amplitude  of  the  traveling  and  stationary  modes  are  the  same.  Owing  to  its  higher 
growth  rate,  traveling  mode  dominates  most  of  the  downstream  development  and  the  growth  of  the 
stationary  mode  is  suppressed.  In  this  case,  energy  cascades  into  (2,2),  (3,3)  etc.,  modes  which  are 
harmonics  of  the  primary  (1,1)  traveling  mode. 

Growth  rates  of  the  stationary  and  traveling  modes  begin  to  depart  from  their  linear  values 
when  the  disturbance  amplitude  reaches  about  4  percent.  As  the  amplitude  increases,  the  primary 
modes  reach  quasi-equilibrium  states.  Large  mean  flow  distortion  caused  by  the  nonlinear 
disturbances  yields  a  skin-friction  value  which  is  significantly  above  the  laminar  value. 

Finally,  we  use  the  two-dimensional  eigenvalue  approach  to  perform  a  secondary  instability 
analysis  of  the  three-dimensional  boundary-layer  flow  modulated  by  the  presence  of  a  nonlinear 
stationary  crossflow  vortex.  We  find  that  this  meanflow  is  subject  to  an  instability  whose  frequency 
is  an  order  of  magnitude  higher  than  the  frequency  of  the  most  amplified  traveling  mode  given  by 
linear  stability  analysis  of  the  boundary-layer  profiles.  A  similar  high  frequency  disturbance  was 
also  observed  in  the  experiments  of  Poll  (1985)  and  Kohama,  Saric  &  Hoos  (1991). 
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Streamwise  (f/,)  and  crossflow  (Uc)  velocity  profiles  in  the  swept  Hiemenz  flow. 


Variation  of  inviscid  streamline  angle  (6)  and  crossflow  Reynolds  number 


Integrated  growth  for  fixed  spanwise  wavenumber  p  = 


0.035 
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Figiire  5.  Growth  rates  (most  amplified)  and  orientation  (y/)  of  the  two  families  of  unstable  modes  at  i?  =  500,  R  =  300. 


Disturbance  growth  rate  for  P=  A  and  R  =  250. 


traveling  disturbance  with  frequency 
Figure  6.  Concluded. 


Disturbance  growth  rate  for  13=  A  and  R  =  500. 


(b)  traveling  disturbance. 

Figure?.  Concluded. 


simulations  of  Streett  (1993)  for  stationary  vortex.  (R  =  500,  /?  =  .4,  growth  rate 


Figure  9.  Computed  wall  vorticity  distribution  in  the  presence  of  nonlinear  stationary  vortices, 
P=A,R  =  500. 
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0.752774 


Variation  of  u  velocity  (for  R  =  500)  in  y-z  plane  in  the  presence  of  stationary  crossflow  vortex.  Two  spanwise  wavelengths  are 
shown.  The  view  is  from  upstream. 


10.  Continued. 


10.  Continued. 


10.  Concluded. 
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Figure  12.  Velocity  vector  plot  at  i?  =  650  projected  onto  a  plane  perpendicular  to  the  vortex 


Variation  of  vortex  axial  velocity  in  the  spanwise  direction  at  i?  =  600. 
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(a)  along  the  vortex 

Figiire  14.  Disturbance  amplitude  functions  for  the  stationary  vortex  (mode  1)  and  its  harmonics,  at  i?  =  600. 
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Figure  14.  Concluded. 
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Figure  15.  Velocity  component  along  the  vortex  at  different  locations  z  across  the  vortex,  R  =  500. 


second  index  refers  to  spanv/ise  wavenumber. 


Figure  16.  Concluded. 


Evolution  of  maximum  (over  y)  amplitude  of  the  disturbance  velocity  components  i.u,v,w)  for  stationary  and  traveling  modes. 
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Variation  of  in  the  y-z  plane  for  the  wave  interaction  case  with  A,  =  .01  percent. 
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19.  Continued. 
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19.  Concluded. 


Fig\ire  20.  Spanwise  variation  of  chordwise  velocity  (M.tationary  and  u^)  aty  =  1.048  for  the  wave  interaction  case  with  A,  =  .01  percent. 


Variation  of  maximum  and  minimum  Ugtationary  and  “rm»  aty  =  1.048  with  Reynolds  number,  u^,  at  locations  where  «,t*tion»ry 
is  maximum  and  minimum  is  also  given. 


Chordwise  iCfu)  and  spanwise  (C^)  skin  friction  coefficients  for  cases  (a),  (b),  (c)  of  figure  16. 


different  numerical  resolutions  are  shown. 
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Figure  25.  The  distribution  of  Xj-total  velocity  component  in  plane  at  -  2.82.  The  lugl  eigenfunction  was  assigned  an  amplitude 
of  5  percent. 
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Abstract 

The  recently  developed  method  of  Parabolized  Stability  Equations  (PSE)  offers  a  fast  and 
efficient  way  of  analyzing  the  spatial  growth  of  linear  and  nonlinear  (convective)  disturbances  in 
shear  layers.  For  incompressible  flows,  the  governing  equations  may  be  represented  either  in 
primitive  variables  or  by  using  other  formulations  obtained  by  eliminating  the  pressure  gradient 
(e.g.,  vorticity/stream-function  formulation).  On  the  other  hand,  for  compressible  flows,  primitive 
variables  offer  a  natural  and  the  only  choice.  We  show  that  the  primitive-variable  fr  ’  vilation  is  not 
well-posed  due  to  the  ellipticity  introduced  by  ^  ^  term  and  marching  solution  aally  blows 
up  for  sufficiently  small  step  size.  However,  it  is  shown  that  this  difficulty  can  be  overcome  if  the 
minimum  step  size  is  greater  than  the  inverse  of  the  real  part  of  the  streamwise  wavenumber,  Or- 
An  alternative  is  to  drop  the  ^/dx  term,  in  which  case  the  residual  ellipticity  is  of  no  consequence 
for  marching  computations  with  much  smaller  step  sizes. 


1.  Introduction 

The  laminar-turbulent  transition  process  in  a  boundary  layer  is  of  great  fundamental  and 
practical  interest  in  fluid  mechanics.  In  many  cases,  transition  occurs  via  the  destabilization  and 
subsequent  growth  of  wave  structures  in  the  boundary  layer.  Classical  theories  concerning  the 
amplification  of  these  waves  use  quasi-parallel  assumption  and  ignore  the  growth  of  the  boundary 
layer.  Other  theories  (Gaster  1974),  including  the  multiple-scales  method  (Saric  &  Nayfeh  1975), 
can  deal  with  the  boundary-layer  growth  locally  for  given  problems.  Direct  Navier-Stokes  solutions 
(Fasel  &  Konzelmann  1990)  could  give  very  satisfactory  results  but  at  the  cost  of  much  more  CPU 
time  and  larger  memory  size.  The  recently  developed  method  of  parabolized  stability  equations 
(PSE)  (Herbert  1991,  Chang  et  al.  1991,  Bertolotti  et  al.  1992)  offers  a  fast  and  efficient  way  to 
analyze  the  spatial  growth  of  instabilities  in  the  boundary  layer.  In  order  to  briefly  describe  the  PSE 
method,  let  us  first  consider  linearized  Navier-Stokes  equations  in  primitive  variable  form: 


du  ..du  dU  .j.du  dU  fjjdu  dU  ^  1 

dt  dx  dbc  ay  dy  dz  dz  dx  R 

d}  dV  ..du  dV  „.d}  dV  dp  1 

--■\.U-—  +  u—+V—  +  v—  +  W—  +  w—+-^  =  — 
dt  dx  dx  dy  dy  dz  dz  dy  R 


d^u  d^u  d^u 
dx^^'dy^^'d^ 


dh  d\ 
dx^  dy 


dz^ 


dw  j.dw  dW  ..  dv  dW  dv  dW  dp  1 

dt  dx  dx  dfy  dy  dz  dz  dz  R 

du  du 
dbc  ^  dy^  dz 


d^w  d^w  d'^w 
dx^ 


(la) 

(lb) 

(10 

(Id) 


where  x,  y  and  z  are  the  streamwise,  wall-normal  and  spanwise  coordinates,  respectively,  and  U,  V, 
W  are  the  corresponding  mean  velocity  components  while  u,  v,  w  are  the  disturbance  velocity 
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components.  Furthermore,  p  is  the  pressure  and  R  is  the  Reynolds  number,  R  =  U^L  /  v ,  where  f7„ , 
L  are  the  reference  scales  for  velocity  and  length  and  v  is  the  kinematic  viscosity. 

In  the  linear  PSE  method,  a  wave-like  disturbance,  0,  representing  velocity,  pressure  or 
vorticity,  etc.,  is  assumed  to  be  of  the  form 

0(jc,  y,z,t)  =  ^x,  y)E  (2a) 


(2b) 


Here  a  and  P  are  the  x  and  ^  wave  numbers,  mis  the  disturbance  frequency  and  $( x,y)  is  the  shape 
function.  The  physical  quantities  {u,v,  etc.)  can  be  obtained  by  adding  the  complex  conjugate 
component.  In  the  above  decomposition  it  has  been  assumed  that  the  mean  flow  is  independent  of 
the  spanwise  coordinate  z.  The  first  and  second  derivatives  of  ip  can  be  written  as 


0,  =(ia0  +  0,)E 

(2c) 

+  2tof0,  +  ^„^E 

(2d) 

II 

(2e) 

(2f) 

Substituting  (2  and  3)  into  the  linearized  Navier-Stokes  equations,  which  could  be  in  primitive 
variable  form,  stream-function/vorticity  form  or  other  forms,  we  obtain  a  set  of  equations  with  the 
$(x,y)'s  and  a  as  unknowns.  Since  there  is  now  one  more  unknown  (namely,  a(x))  than  the 
equations,  another  condition  is  needed  for  the  closure  of  the  system.  We  take  advantage  of  the  slow 
variation  of  the  mean  flow  in  the  streamwise  direction  and  impose  a  condition  on  o(x)  such  that 
"most"  of  the  waviness  and  growth  of  the  disturbance  su’e  absorbed  into  the  exponential  function  E, 
making  the  shape  function  0(x,y)  slowly  varying  with  x.  Hence,  the  term  containing  (in  (3b)) 
can  be  dropped  and  we  surive  at  a  set  of  new  equations  in  which  the  only  second-order  derivatives 
are  those  with  respect  to  y.  We  will  call  these  equations  parabolized  stability  equations  (PSE)  and 
ask  the  question  whether  these  equations  are  indeed  parabolic,  i.e.,  given  appropriate  initial  data 
can  one  find  a  solution  by  marching  edong  the  streamwise  direction? 

When  the  basic  flow  is  two-dimensional,  the  pressure  term  can  be  eliminated  and  the 
governing  equations  can  be  represented  in  stream  function  and  normal  vorticity  formulations.  When 
the  basic  flow  is  three-dimensional  it  becomes  difficult,  particularly  for  nonlinear  disturbances,  to 
reduce  the  equations  to  a  form  other  than  primitive  variable  form.  In  any  case,  for  compressible 
flows,  primitive  variables  offer  a  natural  and  the  only  choice.  Since  our  ultimate  goal  is  to  study  the 
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stability  of  compressible  flow,  we  can  learn  about  the  low  Mach  number  behavior  of  the  compressible 
PSE  equations  by  studying  the  numerical  stability  of  incompressible  PSE  equations  in  primitive 
variable  formulation.  These  equations  can  be  derived  by  using  (1)  and  (2)  and  can  be  written  in  the 
following  form 

Lo^4Li^  +  L2^  +  L30  =  O  (3a) 

where  ^  =  {u,v,w,p)  and  Lq-L^  are  square  matrices  given  in  Appendix  I.  Equation  (3a)  is  to  be 
solved  subject  to  the  constraint  (imposing  the  condition  that  the  shape  function  varies  slowly  with  x  ) 

F(a,^)  =  0  (3b) 

We  will  show  that  when  the  primitive-variables  are  used,  the  equations  are  only  partially 
parabolized.  The  gradient  dp/dt  represents  the  dominant  part  of  the  residual  ellipticity.  For  easy 
reference,  we  call  this  term  as  the  pressure  gradient.  (We  should  stress  here  that  ^/dx  represents 
only  a  small  part  of  the  physical  pressure  gradient  dpidx.  Most  of  the  physical  pressure  gradient  is 
carried  by  the  term  iop,  see  Eq.  (2c)).  This  ellipticity  will  cause  the  marching  procedure  to  fail.  The 
same  difficulties  are  encoimtered  in  the  Parabolized  Navier-Stokes  equations  (PNS)  (Rubin,  1981). 

By  suitable  substitutions,  the  PSE  can  also  be  written  in  the  form  of  a  set  of  first-order  partial 
differential  equations  and  an  additional  functional  F  as 

+<B^  +  I^  =  Q  (4a) 

dx  dy 

F(a,y)  =  0  (4b) 

where  =  ^xy,a)  and  ®  =  T(,xy,a)  are  square  matrices,  I  is  the  identity  matrix  and  yr  is  a  vector.  If 
quasi-parallel  flow  is  assumed,  then  (4)  can  be  reduced  to  the  weU-known  Orr-Sommerfeld  equation. 
Equation  (4a)  is  solved  numerically  by  marching  from  the  initial  station  at  x  =  Xq  with  some  initial 
condition,  y(xo,y).  The  solution  at  x  =  Xo  +  &  is  computed  with  a(xo +&)  =  a(xo)  as  a  first 
approximation,  then  a  new  a  is  calculated  using  Eq.  (4b).  Equation  (4a)  is  solved  again  with  the  new 
value  of  oc.  This  process  continues  until  the  solution  converges.  The  marching  is  then  carried  to  the 
next  x-station.  The  stability  of  the  marching  procedure  depends  on  the  discretization  scheme,  the 
iterative  process  for  evaluating  a  and,  most  importantly,  the  mathematical  nature  of  PSE. 

In  the  primitive-variable  form,  the  matrices  in  Eq.  (4a)  for  a  flat  geometry  are  given  as  follows 
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Ji  = 
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0 

0 

0 

0 

-1 
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ia 
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ip 

0 

0 
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0 

0 

0 

0 

-1 

-iaV  +  Vx 

d  +  Vy 

-ipv 

0 

ia 

R 

R 

UyR 

0 

-ioR 

-VR 

0 

-W^R 

-WyR 

-AR 

-iPR 

0 

-VR 

p^-i 

da 

-io)  +  icdJ  +  ipw  + - 

R 

dx 

0  0  0 

0 

0 

O' 

10  0 

0 

0 

0 

0  0  0 

0 

0 

0 

-V  r  0 

0 

1 

R 

0 

-moo 

-R 
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0 

(5) 


(6) 


2tcc 

where  F  =  U  — — ,  U,  V  and  W  are  the  three  velocity  components  of  the  basic  flow,  a  and  P  are  the 
R 

wave  numbers  in  the  x  and  y  -  directions  respectively,  and 


=  {u,v,w,p,dul  dy,du} !  dyf 


(7) 


The  functional  F{a,  v^)  =  0 ,  can  be  chosen  in  several  ways.  For  example,  we  can  impose  the 
condition  that  the  maximum  of  the  velocity  component,  u ,  is  constant  or  that 


_o _ 

•a 


(8) 


where  q  =  [u,v,w)  and  *  represents  the  complex  conjugate.  An  iterative  procedure  for  a  based  on  Eq. 
(8)  is  given  as  follows. 
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(9) 


9 _ 

•0 

0 


We  note  that,  since  a  is  complex,  Eq.  (9)  represents  a  two-dimensional  iterative  map. 

2.  Mathematical  Nature  of  PSE 
For  an  initial  value  problem 

du  .dUn 

uM  =  ny) 


(10) 

(11) 


where  f,ueC'‘,  A,BeC"*'‘  and  ye 91,  the  solutions  can  be  obtained  by  Fourier  transform.  In 
Fourier  space,  Eq.  (10)  and  (11)  become 


^0)  =  /(n) 


(12) 

(13) 


where  u  and  f  are  Fourier  transforms  of  u  and  /,  respectively,  and  is  the  wavenumber  in  y- 
direction.  The  solution  to  (12)  and  (13)  is  of  the  form 


u-g(v)e^ 


(14) 


where  g(T})eC''-  Substitution  of  (14)  into  (12)  yields  an  algebraic  equation  for  the  eigenvalue  A  = 

|t7jA  +  jB-Ar|  =  0  (15) 

Assuming  that  aU  n  eigenvalues  of  Eq.  (15)  have  distinct  corresponding  eigenvectors,  gj(rj),  then 

ft 

u  =  ^AjgjCti)e^J* 
j=i 

The  constants  Aj  e  91  can  be  determined  by  solving  the  nxn  linear  system 

^gj(V)Aj  =  /(??) 

J=i 

Hence,  we  obtain  the  solution  to  Eq.  (10)  and  (11) 
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(16) 


Aje‘’>^gj{v)e^^drt  . 


If  is  unbounded  for  (rjj  ^  oo,  the  integral  in  Eq.  (16)  may  not  exist,  and  the  initial  value  problem 
defined  by  Eqs.  (10)  and  (11)  are  said  to  be  ill-posed  (Kreiss  &  Lorenz  (1989)). 

In  finding  the  eigenvalues  both  matrices  A  and  B  are  used  in  Eq.  (16).  In  fact,  subject  to 
the  condition  that  the  eigenvalues  of  A  are  non-zero  and  distinct,  the  well-posedness  of  Eq.  (10)  is 
determined  by  the  eigenvalues  of  A  only.  This  can  be  easily  verified.  We  decompose  A  such  that 
A  =  PAP~^,  where  P  is  a  matrix  whose  columns  consist  of  the  eigenvectors  of  A,  and  A  is  a  diagonal 
matrix  whose  elements  are  eigenvalues  of  A.  Multiplying  Eq.  (15)  by  [P"^|P1,  we  obtain 

|t7jA  +  P‘^BP-;i/|  =  0 

As  the  matrix  ii]A  +  P~^BP  becomes  diagonally-dominant.  Applying  the  Gershgorin  Circle 

Theorem  (see  Golub  &  Van  Loan,  1983),  we  find  that  the  eigenvalue,  kj ,  approaches  ijOj ,  where  Oj 
is  the  jth  eigenvalue  of  A.  Therefore,  the  eigenvalues  of  A  determine  the  well-posedness  of  Eq.  (10). 
This  is  reminiscent  of  the  classification  of  partial  differential  equations  with  one  dependent  variable 
where  only  the  "principal  part"  is  required.  However,  we  must  stress  that,  if  some  of  the  eigenvalues 
of  A  are  zero,  then  aU  matrices  in  Eq.  (15)  must  be  used  to  obtain  the  correct  result. 

For  various  PSE  formulations,  the  equations  cannot  always  be  written  in  the  form  of  Eq.  (10), 
e.g.,  when  matrix  CB  in  Eq.  (4)  is  singular.  Here,  we  adopt  the  approach  whereby  Eq.  (14)  is 
substituted  into  the  given  form  of  PSE  and  all  matrices  are  kept  for  determining  the  well-posedness. 

We  now  consider  and  analyze  the  well-posedness  of  the  initial  value  problem  associated  with 
the  PSE  equations.  The  system  consisting  of  Eqs.  (3a)  and  (3b)  (or  the  corresponding  system  4(a) 
and  4(b))  is  nonlinear  since  a  appears  in  the  coefficient  matrices.  In  order  to  simplify  the  numerical 
stability  analysis  of  the  marching  procedure,  we  will  assume  that  a  is  known  a  priori.  We  choose  a 
simple  two-dimensional  basic  flow  (W  =  0,  d!  dz  =0)  with  constant  velocity-components  and  apply 
Fourier  transform  in  y-direction,  i.e.,  let  uix,y)  =  u(x)exp(irjy).  After  some  algebra,  we  can  write  Eq. 
(3a)  as  follows. 


d_ 

dx 


-ia 

-IT} 

0  ■ 

V 

0 

C 

V 

-C  +  iaD 

D 

iriD 

D 

-ia 

(17) 


doc  ^  ^ 

where  C  =  iaf/  -ico-i  +  itjV  +  and  D  =  U -2i^.  The  three  eigenvalues  of  the  matrix  in 

dx  R  R  R 

Eq.  (17)  are  given  by 


ki  =  T]-ia,  k^=-r]-ia, 
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(18) 


For  a  well-posed  initial  value  problem,  it  is  required  that  Real  Ay  <  0  for  |7jj  ->  «> .  Since  jj  can  take 
any  value  for  general  initial  conditions,  Ai  and  A2  have  positive  real  parts  for  some  tj.  Therefore,  Eq. 
(17)  is  not  well-posed  for  initial  value  problem  and  some  approximation  must  be  made  for  a  stable 
marching  solution.  We  will  discuss  it  in  the  next  section. 

Due  to  the  complicated  nature  of  A3,  it  is  not  easy  to  find  the  sign  of  its  real  part.  However,  in 
order  to  determine  the  well-posedness,  we  only  need  to  know  the  behavior  of  A3  as  lij]  ,  provided 
that  the  coefficient  of  the  highest  order  term  in  rj  is  not  purely  imaginary.  In  this  limit,  we  have 


A3  — »  — ■ 


R  {U-2iaJR) 


Therefore, 


{U  +  2ailR)  +  2ia,IR 
R  {U  +  2ai/Rf+4af/R^ 


(19) 


R6fll^A3^  — >  “ 


_ U  ^lajlR _ 

^  {U  Rf  +4al !  R"^ 


(20) 


Hence,  for  growing  waves  (o^  <  0),  real  A3  >  0  for  sufficiently  small  17,  i.e., 

[7  <^5-  (21) 

For  large  Reynolds  number  and  small  0%  (typically  0(10^)  and  0(10-*),  respectively),  condition  (21)  is 
satisfied  only  in  a  very  small  neighborhood  of  the  wall.  As  we  will  show  later,  difficulties  in 
numerical  solution  never  arise  from  A3. 


3.  Stability  of  the  Marching  Procedure 

We  will  now  show  that  marching  procedure  can  be  made  stable  by  using  a  sufficiently  large 
step  size,  &,  as  in  PNS  for  mean  flow,  and  we  will  derive  a  condition  for  stability.  We  note  that, 
since  the  eigenvalues  are  distinct,  the  matrix  in  Eq.  (12)  is  diagonalizable,  and  therefore,  Eq.  (17)  can 
be  written  as 


0  0 
A2  0 
0  A3 


(22) 


where  ^  is  a  vector  whose  components  are  linear  combinations  of  u,  v  and  p .  These  are  three 
decoupled  equations,  each  of  which  is  of  the  form 


d(p 

dx 


=  A^ 


(23) 


If  this  equation  is  solved  with  backward  difference  with  step-size  &c,  von  Neumann  analysis  leads  to 
an  amplification  factor  given  by 


(24) 


1 

1-SxX 


Numerical  stability  requires  |y(<l,  which  can  be  satisfied  if  (l-X^Sxf  >1  and  hence, 

dx  >  2X^  +  ^)-  Since  X  depends  on  rj,  then  minimum  step  size  is  given  by 


(25) 


Equation  (25)  implies  that  there  is  no  step-size  restriction  if  A,  <  0,  otherwise  the  step  size  has  to  be 
greater  than  some  finite  value  if  the  marching  is  to  be  stable.  For  example,  von  Neumann  analysis 
of,  say,  the  equation  associated  with  Xq  leads  to  an  amplification  factor  given  by 


^2  = 


^1-^Tf-a.  I&j  +  ia^Sx 


(26) 


Numerical  stability  requires  |y2|<  1  for  all  tj-  This  can  be  achieved  only  if  &  >  j-i-r.  For  smaller  Sx, 

r*r| 

numerical  instability  occurs  for 


V 


This  is  graphically  depicted  in  Fig.  1  for  real  a 

If  Eq.  (4a)  is  discretized  directly  by  using  backward  difference  for  x-derivatives  and  central 
difference  for  the  y-derivatives,  then  Von  Neumann  analysis  of  numerical  stability  leads  to  the 
following 
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where  y  is  the  amplification  factor, 
^  =  1+  iaX .  The  six  roots  of  Eq.  (27)  are 
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(30) 


,  ,  -ia±^f-a^ 

The  amplification  factor  associated  with  is  given  by 
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76=- 
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„2  sin  *(7/^) 


Sy^ 


+  ia^& 


which,  again,  leads  to 
A-  -JL 

■  K1 


(31) 


(32) 


for  numerical  stability.  The  minimum  step  size  as  given  by  Eq.  (32)  implies  that  a  maximum  of  2;r 
steps  per  disturbance  wavelength  are  allowed  for  the  marching.  This  is  consistent  with  the  behavior 
of  the  numerical  solution  first  reported  by  Chang  et  al.  (1991).  We  should  note  that,  in  PSE,  we  solve 
for  the  slow  varying  shape  functions  and  so  this  step-size  restriction  does  not  cause  problems  in 
terms  of  accuracy;  indeed,  the  numerical  examples  given  in  Joslin  et  al.  (1993)  show  that  the  PSE 
solution  obtained  with  3  steps  per  wavelength  for  ToUmien-Schlichting  waves  in  a  Blasius  boundary 
layer  is  in  excellent  agreement  with  very  accurate  direct  simulation  of  Navier-Stokes  equations  using 
60  grid  points  per  TS  wavelength. 

When  step-size  smaller  than  =  l/|arl  is  required  for  either  higher  resolution  or  for  the 
convergence  of  nonlinear  terms  in  an  implicit  numerical  scheme,  a  further  approximation  is 
sometimes  made,  i.e.,  dropping  the  pressure  gradient  term  ^/dx  (see  Chang  et  al.  1991).  We  point 
out,  however,  that  most  of  the  physical  pressure  gradient  is  absorbed  into  the  term  iqp,  and  that 
c^/dc  is  very  small  in  comparison  (refer  to  Eq.  2(c)).  After  ^/obc  is  dropped,  Eq.  (17)  can  be  reduced 
to  a  set  of  2  equations  by  eliminating  p .  Canying  out  the  well-posedness  analysis,  we  obtain  two 
eigenvalues 


Aj  = 


A2  -  - 


C 

D 


where  C  and  D  are  the  same  as  those  in  Eq.  (18).  Here,  gives  rise  to  ill-posedness  due  to 
viscosity,  which  is  insignificant  as  explained  before.  We  note  that  represents  ill-posedness  if  a;  < 
0  (i.e.,  for  growing  waves).  Using  Eq.  (25),  we  obtain  the  minimum  step-size  required  for  a  stable 
marching  scheme. 
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In  a  typical  problem,  a,  is  an  order  of  magnitude  smaller  than  ■  Therefore, 


For  a  problem  in  which  /|al  is  small,  we  see  that  the  restriction  on  step-size  is  relaxed  by  at  least 
an  order  of  magnitude  in  comparison  with  Eq.  (32). 


4.  Application  to  Blasius  Boundary  Layer 

The  above  analysis  was  simplified  in  that  the  mean  flow  (U,V)  was  considered  constant.  We 
now  apply  the  results  of  the  above  analysis  to  a  realistic  problem,  i.e.,  the  boundary-layer  flow  on  a 
flat-plate.  The  mean  flow  in  this  case  is  governed  by  the  Blasius  equation 


^  1  ^ 

de  r  de 


=  0 


/■(0)  =  r(0)  =  0;  r(4)-»l  as 


u=r 


where  v  is  the  kinematic  viscosity,  f/,  the  boundary-layer  edge  velocity  and  ^  the  similarity  variable 


and  Xq  and  yo  baing  the  dimensional  values  ofx  andy. 

The  Blasius  boundary-layer  is  subject  to  ToUmien-Schlichting  (TS)  instability.  We  consider  a 
disturbance  of  fixed  frequency  F  =  0.7  x  10'*  where 


f  being  the  dimensional  frequency.  For  this  disturbance,  the  initial  wavenumber  at  i?  =  500  is  = 
.106.  As  the  wave  travels  downstream,  the  wavenumber  gradually  changes  to  approximately  = 
.103. 
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Calculations  were  started  at  i?  =  U^LJv  =  500,  where  L  =  ^vxqIU^  .  According  to  Eq.  (32),  the 
minimum  step-size  is  =9.5  based  upon  the  length  scale  at  R  =  500.  The  fourth-order  accurate 
compact  difference  scheme  of  Malik  et  al.  (1982)  was  used  for  wall  normal  discretization  and  one  step 
backward  Euler  discretization  was  employed  in  the  streamwise  direction.  We  ran  our  PSE  code  with 
term  included  and  Ax  =  11  (greater  than  Ax^„).  The  growth  rate  plotted  in  Figure  2(a)  shows 
a  smooth  solution.  The  growth  rate  in  this  figure  has  been  made  nondimensional  with  local  length 
scale  which  varies  with  x;  however,  a  constant  length  scale  was  used  in  the  marching  computations. 
Next  we  use  Ax  =  8  (smaller  than  Ax^i^),  the  solution  blows  up  shortly  after  the  marching  starts. 
Now,  we  drop  $/dx  and  further  reduce  the  step-size  to  5,  a  smooth  solution  is  again  obtained. 
These  calculations  confirm  the  estimate  for  the  minimum  step-size  and  show  that  the  dropping  of 
^/dx  relaxed  the  step-size  restriction.  We  present  additional  calculations  in  Figure  2(b)  where  Ax 
=10  is  used  with  ^/dx*Q  and  the  solution  is  smooth.  However,  with  Ax  =  9,  the  solution  blows  up 
towards  the  end  of  the  computation.  This  shows  that  the  stability  condition  derived  in  the  previous 
section  provides  a  remarkably  good  prediction  of  the  behavior  of  the  PSE  solution.  With  dp! dx.  =  Q, 
stable  solution  can  be  obtained  for  smaller  Ax.  The  comparison  of  the  solution  with  Ax  =  5  and  10 
shows  that  the  approximation  ^tdx  =  0  does  not  introduce  any  error,  at  least  in  this  case. 

Figure  3  shows  the  growth  rate  results  for  the  same  frequency  disturbance  and  ^/A:  =  0. 
Three  different  step  sizes  were  chosen:  Ax  =  .3,  .25  and  .2.  With  Ax  =  .3,  the  solution  is  generally 
smooth  but  appears  to  develop  some  very  slight  wiggles  near  the  end.  However,  with  Ax  =  .25  and  .2 
solution  blows  up.  The  numerical  instability  occiu's  earlier  for  the  smaller  step  size. 

The  Ax  =  .2  calculation  is  repeated  with  ia/R  dropped  to  eliminate  elliptic  effect  due  to 
viscosity.  The  behavior  of  the  solution  is  unaffected.  This  indicates  that  these  oscillations  are  not 
caused  by  ellipticity  due  to  the  viscous  terms.  Finally,  we  eliminate  the  iteration  for  a,  and  do  the 
calcidation  with  fixed  real  a  (the  value  at  /{  =  500).  The  growth  rate  shown  in  Figure  3  is  inaccurate, 
but  the  growth  rate  curve  is  smooth.  This  suggests  that  the  oscillations  may  have  been  caused  by 
the  nonlinear  iterative  process  for  determination  of  a.  However,  one  cannot  be  absolutely  sure  in 
view  of  the  stability  condition  (33).  For  the  present  case  of  =  .0055  and  ®  •!>  i*-  follows 

fi'om  (33)  that  Ax^j,,  =  .13  which  is  not  far  from  the  value  of  .25  or  .2  used  in  Figs.  3(b)  and  (c), 
respectively.  The  fact  that  smooth  solution  was  obtained  in  Fig.  3(c)  with  a  fixed  does  not 
necessarily  suggest  that  (33)  is  not  operative  since  retd  (a)  was  used  and  with  =  0  Eq.  (33)  gives 

d*niin  ~  0  • 

In  order  to  clarify  whether  (33)  is  operative,  we  perform  additional  computations  and  report 
the  results  in  Fig.  3(d)  which  contains  the  correct  solution  from  Fig.  2(a)  as  well  as  the  smooth 
solution  from  Fig.  3(c)  with  a,-  =  0  and  Ax  =  .2.  Also  included  are  solutions  using  fixed  a  =  .104  -  i 
.001  and  a  =  .104  -  i  .004  ( p,  and  iodR  have  been  dropped).  Equation  (33)  for  =  0,  =  -.001  and 

-.004  yield  minimum  step  size  of 0,  .023  and  .092,  respectively.  In  the  first  two  cases,  the  step  size  of 
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Ax  =  .2  used  in  the  computation  is  much  larger  than  the  critical  value  given  by  (33)  and  therefore  the 
solution  is  smooth.  Since  initially  the  value  of  a,  is  close  to  the  correct  solution  given  by  (4b),  the 
computed  growth  rate  is  close  to  the  correct  growth  rate.  However,  further  downstream  the 
computed  growth  rate  is  in  error  since  (4b)  is  not  satisfied.  With  =  -.004,  the  computational  step 
size  is  not  very  far  from  the  critical  value  and,  therefore,  solution  becomes  oscillatory  and  eventually 
blows  up.  In  this  case,  since  initially  «£  =  -.004  is  not  close  to  the  solution  of  (4b)  the  growth  rate  is 
in  error  at  small  Reynolds  numbers.  However,  as  the  correct  growth  rate  approaches  the  value  of 
-004,  the  solution  with  Ac  =  .2  and  without  imposition  of  4(b)  approaches  the  correct  solution  even 
though  it  is  oscillatory  due  to  the  numerical  instability.  This  example  clearly  suggests  that 
numerical  instability  observed  with  Ax  =  .2  is  associated  with  the  condition  given  by  (33). 

We  point  out  that  the  values  of  Ax  =  .3  and  .2  correspond  to  about  200  to  300  marching  steps 
per  TS  wavelength.  Clearly,  in  no  linear  or  nonlinear  applications  one  needs  to  take  such  a  small 
step  and,  therefore,  Eq.  (33)  is  of  no  significance  for  most  practical  computations. 

5.  Conclusions 

The  mathematical  nature  of  parabolized  stability  equations  (PSE)  is  studied.  It  is  shown  that 
the  primitive  variable  formulation  is  mathematically  ilbposed  due  to  the  pressure-gradient  term. 
The  condition  for  the  stable  marching  solution  is  derived.  Examples  from  the  linear  stability  of  two- 
dimensional  Blasius  boundary-layer  are  given  to  show  that  this  condition  gives  reasonable  estimates 
of  the  numerical  behavior  of  the  parabolized  equations.  The  results  also  show  that  the  ellipticity 
associated  with  the  viscous  term  is  insignificant. 
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Appendix  I 
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Effect  of  marching  step  size  on  the  growth  rate  of  the  TS  disturbance  of  frequency  F  =  .7  x  10"^ .  The  critical  step  size  for  stable 
solution  is  =  9.5. 


3(c).  Continued. 


